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Abstract 

We  apply  an  inverse  problem  formulation  to  determine  characteristics  of  a  defect  from  a  perturbed 
electromagnetic  interrogating  signal.  A  defect  (crack)  inside  of  a  dielectric  material  causes  a  disruption, 
from  reflections  and  refractions  off  of  the  interfaces,  of  the  windowed  interrogating  signal.  We  model  the 
electromagnetic  waves  inside  the  material  with  Maxwell’s  equations.  Using  simulations  as  forward  solves, 
our  Newton-based,  iterative  optimization  scheme  resolves  the  dimensions  and  location  of  the  defect. 
Numerical  results  are  given  in  tables  and  plots,  standard  errors  are  calculated,  and  computational  issues 
are  addressed. 


1  Introduction 

The  problem  we  are  considering  is  detecting  a  crack  inside  of  a  dielectric  material  using  electromagnetic 
interrogation.  The  idea  is  to  measure  the  reflected  and/or  transmitted  signals  and  use  the  data  to  solve  an 
inverse  problem  to  determine  certain  characteristics  of  the  crack,  e.g.,  location  and/or  width.  Possible  appli¬ 
cations  of  this  procedure  include  quality  assurance  in  fabrication  of  certain  dielectric  materials,  or  damage 
detection  in  aging  materials  for  safety  concerns.  Further  applications  of  electromagnetic  interrogation,  as 
well  as  additional  problem  formulations  and  solutions,  can  be  found  in  [2]. 

First  we  simplify  the  problem  by  considering  a  linearly  polarized,  pulsed  interrogating  signal  which  reduces 
the  problem  to  one  spatial  dimension.  We  start  by  defining  an  inverse  problem  for  determining  the  crack’s 
dimensions.  We  assume  that  we  have  data  from  sensors,  located  in  front  of  and/or  behind  the  material,  that 
receive  the  electromagnetic  signal  after  it  is  reflected  from  (or  passes  through)  the  material.  We  then  compute 
simulated  signals  with  approximations  to  the  crack’s  characteristics  and  apply  an  optimization  routine  to  a 
modified  least  squares  error  between  this  simulated  signal  and  the  given  data.  In  our  computations  we  use 
Maxwell’s  equations  and  a  Debye  polarization  equation  to  model  the  signal,  and  solve  these  equations  using 
a  Finite  Element  method  in  space  and  Finite  Difference  methods  in  time.  Thus  the  optimization  routine 
finds  those  crack  characteristics  which  generate  a  simulated  signal  that  most  closely  matches  the  given  data. 
In  this  sense  we  determined  an  estimate  to  the  “true”  crack  characteristics. 

In  Section  2  we  define  the  equations  that  we  have  chosen  in  order  to  model  the  electromagnetic  waves 
inside  the  material.  We  further  distinguish  between  two  distinct  problem  types  that  we  will  address,  ( Problem 
1  and  Problem  2).  Section  3  contains  the  details  of  our  numerical  methods  for  the  simulations.  We  introduce 
the  inverse  problem  formulation  for  Problem  1  in  Section  4,  and  later  improve  upon  it  in  Section  4.2. 
Numerical  results  of  the  inverse  problem  are  displayed  in  Section  4.3. 

In  Section  5  we  begin  addressing  Problem  2.  Similarities  and  differences  between  the  computational 
issues  between  the  two  problems  are  pointed  out.  A  more  sophisticated  optimization  method  is  described 
in  Section  5.3  and  associated  numerical  results  are  given  in  Section  5.5.  Sections  5.5.1  and  5.5.2  explore  the 
effects  of  adding  random  noise  to  the  data,  both  relative  and  constant  variance.  In  the  latter,  we  compute 
standard  error  estimates. 


2  Problem  Description 

We  interrogate  an  (infinitely  long)  slab  of  homogeneous  nonmagnetic  material  by  a  polarized,  windowed 
signal  (see  [2]  for  details)  in  the  Thz  frequency  range  (see  Figure  1).  We  start  with  a  wave  normally  incident 

1  Center  for  Research  and  Scientific  Computation,  Raleigh,  NC 
2NASA  Langley,  VA. 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
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Figure  1: 


Problem  1 :  Dielectric  slab  with  a  material  gap  in  the  interior.  Possible  sensors  in  front  and  behind. 


Figure  2:  The  domain  of  the  material  slab:  O  =  [z i,  24]. 


on  a  slab  which  is  located  in  fl  =  [zi,  z4],  (  0  <  Zi  <  Z4  <  1)  with  faces  parallel  to  the  x-y  plane  (see  Figure 
2).  We  denote  the  vacuum  outside  of  the  material  to  be  flo-  The  electric  field  is  polarized  to  have  oscillations 
in  the  x-z  plane  only.  Restricting  the  problem  to  one  dimension,  we  can  write  the  electric  and  magnetic 
fields,  E  and  H  respectively,  as  follows 

E{t,  x)  =  iE(t,  z) 

H(t,x)  =  jH(t,  z), 

so  that  we  are  only  concerned  with  the  scalar  values  E(t,  z)  and  H(t ,  z). 
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Maxwell’s  equations  become: 

(la) 

(lb) 

where  D(t,  z)  is  the  electric  flux  density,  /to  is  the  magnetic  permeability  of  free  space,  a  is  the  conductivity, 
and  Js  is  a  source  current  density  (determined  by  our  interrogating  signal).  We  take  the  partial  derivative 
of  Equation  (la)  with  respect  to  z,  and  the  partial  of  Equation  (lb)  with  respect  to  t.  Equating  the 
terms  in  each,  and  thus  eliminating  the  magnetic  field  H ,  we  have: 

E"  =  /to  (jd  +  ctE  +  js ^  , 

(where  '  denotes  2  derivatives  and  '  denotes  time  derivatives) . 

Note  that  we  have  neglected  magnetic  effects  and  we  have  let  the  total  current  density  be  J  =  Jc  +  Js, 
where  Jc  =  aE  is  the  conduction  current  density  given  by  Ohm’s  law. 

For  our  source  current,  Js,  we  want  to  simulate  a  windowed  pulse,  i.e.,  a  pulse  that  is  allowed  to  oscillate 
for  one  full  period  and  then  is  truncated.  Further,  we  want  the  pulse  to  originate  only  at  2  =  0,  simulating 
an  infinite  antenna  at  this  location.  Thus  we  define 


dE 


dH 


dz  dt 


dH 

dz 


dD 

~dt 


+  it  E  +  Js 


Js(t,z)  =  S(z)sin(u;t)I[0jtf](t ) 

where  u>  is  the  frequency  of  the  pulse,  tf  =  27r/o>  is  fixed,  /[o,t,](i)  represents  an  indicator  function  which  is 
1  when  0  <  t  <  tf  and  zero  otherwise,  and  d(z)  is  the  Dirac  delta  distribution. 

Remark  1  Computationally,  having  a  truncated  signal  introduces  discontinuities  in  the  first  derivatives 
which  are  not  only  problematic  in  the  numerical  simulations  (producing  spurious  oscillations),  but  are  also 
essentially  non-physical.  Therefore  in  our  implementation  we  actually  multiply  the  sine  function  by  an 
exponential  function  (see  Figure  3)  rather  than  the  traditional  indicator  function.  The  exponential  is  of  the 
form 


where  6  controls  the  horizontal  shift,  a  determines  the  width  and  (3  determines  the  steepness  of  the  sides. 
A  value  of  9  =  2  provides  a  sufficient  buffer  before  the  signal  to  avoid  leading  oscillations  due  to  the  initial 
discontinuity.  For  notational  consistency  we  will  continue  to  denote  this  function  as 

The  electric  flux  density  inside  the  material,  given  by  D  =  eoCoo E  +  P,  is  dependent  on  the  polarization, 
P.  Note  that  e0  is  the  permittivity  of  free  space  and  is  the  relative  permittivity  in  the  limit  of  high 
frequencies  (thus  the  notation  of  00).  For  computational  testing  we  assume  for  this  presentation  that  the 
media  is  Debye  and  thus  we  use  the  following  polarization  model  inside  O: 

TP  +  P  =  e0(es  -  Coo )E 

where  es  is  the  static  relative  permittivity  and  r  is  the  relaxation  time.  We  also  assume  P(0,  z)  =  0.  Note 
that  in  the  vacuum  outside  of  f2,  P  =  0. 

In  order  to  represent  D  in  the  entire  domain,  we  use  the  indicator  function  Iq  which  is  1  inside  12  and 
zero  otherwise.  Thus 

D  =  eo  E  +  eo(eoo  —  1  )Iq.E  +  IqP- 

In  order  to  have  a  finite  computational  domain,  we  impose  absorbing  boundary  conditions  at  2  =  0  and 
2  =  1,  which  are  modelled  as 


E-cE' 
E  +  cE' 


=  0 


J  z= 0 


=  0. 
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Figure  3:  Our  choice  of  a  smooth  indicator  function  and  the  resulting  (smoothly  truncated)  interrogating 
signal.  In  this  example  9  =  2,  a  =  .9,  and  f3  =  10. 


With  these  boundary  conditions,  any  incident  signal  passes  out  of  the  computational  domain,  and  does  not 
return.  We  model  this  by  forcing  it  to  be  absorbed  by  the  boundary.  Also  we  assume  zero  initial  conditions, 
i.e., 


E(0,z)  =  0 
E(0,z)  =  0. 

Thus  our  entire  system  can  be  written 

Moeo(l  +  (e<x>  ~  1  )Iq)E  +  no IqP  +  Ho&E  —  E"  =  —fio Js 

TP  +  P  =  e0(es  -  )E 
[E  -  cE']z= o  =  0 
[E  +  cE']z=1  =  0 
E(0,z)  =  0 
E(0,z)  =  0 


in  17  U  Oo 
in  fl 


with 

Js(t,z)  =  6(z)sin(ujt)I[o,tf](t). 

In  this  formulation  we  have  initially  assumed  a  single  slab  of  a  dielectric  contained  in  O  =  [21,24].  Thus 
Iq  —  1  if  21  <  2  <  24,  and  zero  otherwise.  We  now  introduce  a  gap  (crack)  consisting  of  a  vacuum  in  the 
interior  of  the  material  as  depicted  in  Figure  4.  If  the  gap  is  (22,23)  then  we  redefine  O  =  {2I21  <  2  < 
22  or  Z3  <  2  <  24}.  We  will  refer  to  this  formulation  as  Problem  1  (recall  Figure  1).  Later  we  will  discuss 
a  second  formulation,  Problem  2,  where  the  gap  is  between  the  dielectric  slab  and  a  metallic  (conducting) 
backing,  as  shown  in  Figure  5.  This  will  require  slightly  different  boundary  conditions  (reflecting  instead  of 
absorbing  at  2  =  1,  where  the  metal  backing  begins),  but  otherwise  the  numerical  methods  and  analysis  are 
the  same. 
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Figure  4:  The  domain  of  the  material  slab  with  an  interior  gap  between  z 2  and  23:  O  =  {z\zi  <  z  < 

Z2  or  23  <  z  <  24}  . 


x 


Figure  5:  Problem  2:  Dielectric  slab  and  metallic  backing  with  a  gap  in  between.  Possible  sensors  only  in 
front. 


3  Numerical  Solution 

3.1  Finite  Elements 

We  apply  a  Finite  Element  method  using  standard  linear  one  dimensional  basis  elements  to  discretize  the 
model  in  space.  Let  N  be  the  number  of  intervals  in  the  discretization  of  z,  and  h  =  1/N,  then  the  Finite 
Element  discretization  has  an  order  of  accuracy  of  0(h2).  For  implementation  we  scale  time  by  t  =  ct  and 
the  polarization  by  P  =  P/cq  for  convenience.  The  resulting  system  of  ordinary  differential  equations  after 
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the  spatial  discretization  is  the  semi-discrete  form 


erMe  +  Mnp  +  (770  <tM  +  D  +  B)e  +  Ke  =  ijoJ  (2a) 

Mnp  +  XMnp  =  edXMne,  (2b) 


where  er  =  (1  +  (e^  —  l)/n),  ed  =  es  —  £001  X  =  1  ,  and  770  =  yVo  Ao  ■  Also  e  and  p  are  vectors  representing 
the  values  of  E  and  P  respectively  at  the  nodes  z,  =  ih.  The  mass  matrix  M  has  entries 


Mi:j  =< 


>■=  4>i4>jdz 


where  are  the  basis  functions  (Mn  is  the  mass  matrix  integrated  only  over  fi),  while  the  stability 

matrix  K  has  entries 

Kij  =<  $,</>'•  >■=  fii&'jdz. 

Jo 

The  matrices  D  and  B  result  from  the  boundary  conditions  where 


£1,1  =  1 

Bn+i,n+i  =  1 

and  all  other  entries  are  zero.  Finally,  J  is  defined  as 

Ji  —  0  i  3  J  s  P  •  —  / 


Js(j>idz. 


Note  that  by  differentiating  (2b)  we  can  substitute  into  (2a)  and  obtain  an  equation  only  dependent 
explicitly  on  P  (two  substitutions  are  required  to  eliminate  P  and  P ): 

erMe  +  (rjocrM  +  D  +  B  +  edX  M)e  +  ( K  —  edX  2M)e  +  X2Mnp  =  rj0J. 

Using  shorthand  we  can  write  our  entire  coupled  system  as 

M 1  e  +  M2e  +  M^e  +  X  2p  =  po  J 


p+  Xp=  edXMne 


(3a) 

(3b) 


where  p  =  M(‘lp.  It  is  important  to  mention  that  each  matrix  is  tridiagonal  due  to  the  choice  of  the  linear 
finite  elements. 


3.2  Finite  Differences 

In  order  to  solve  the  semi-discrete  form  of  our  equations  we  consider  two  distinct  finite  difference  methods. 
In  the  first  method  we  convert  the  coupled  second  order  system  of  equations  into  one  larger  first  order  system 
and  simply  apply  a  theta  method  (unless  otherwise  stated,  we  use  6  =  |).  In  the  second  method  we  solve 
first  for  the  polarization  with  a  forward  differencing  scheme  using  the  initial  conditions  and  then  use  that 
to  update  a  second  order  central  difference  scheme  for  the  magnitude  of  the  electromagnetic  field.  We  then 
continue  this  process  iteratively,  alternating  between  solving  for  P  and  for  E. 

Both  methods  are  second  order  in  time  and  space  for  appropriately  smooth  data  (and  with  At  =  O(h)). 
We  have  compared  the  errors  and  the  run  times  of  the  two  methods  for  several  smooth  test  problems  and 
have  determined  that  the  second  method  is  as  accurate,  but  twice  as  fast,  as  the  first  in  all  cases  primarily 
due  to  the  fact  that  the  linear  system  is  of  smaller  dimension  in  the  second  method.  Also,  the  first  method 
incidentally  solves  for  e  in  addition  to  e  and  p,  which  is  superfluous. 

In  our  second  method  we  use  a  second  order  central  difference  scheme  to  solve  (3a).  Thus  we  must  first 
find  an  approximation  to  E(ti,z)  where  t,;  =  iAt.  Note  that  approximating  E  with  its  Taylor  expansion 
around  =  0  and  applying  the  initial  conditions  and  ODE,  one  obtains 

At2 

E(ti,z)  «  -—fj,oJa(0,z). 
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Our  approach  is  to  first  solve  for  p  using  a  0-method,  and  then  use  that  approximation  to  solve  for  e  at 
the  next  time  step.  Thus,  our  finite  difference  approximation  for  (3b)  is 

Vn+i  =Pn+  :  +  ^0  ( edMnen+e  -  pn)  (4) 

where  [en]j  =  E(tn,Zj),  [pn\j  =  MnP(tn,  Zj),  Zj  =  jh,  and  en+g  =  0en  +  (1  —  0)en+ 1  is  a  weighted  average 
of  en  and  en+i  for  relaxation  to  improve  the  stability  of  the  method.  Once  we  have  pn+\  we  can  solve  for 
e„+2-  Applying  second  order  central  difference  with  averaging  to  (3a)  gives 

M2^n+2  =  Men  +  M\en+\  +  At2r]oJn+i  ~  A“At“pn+i.  (5) 

Note  that  in  this  case  M2  is  tridiagonal  and  the  matrix  is  the  same  for  each  time  step,  so  we  may  store 
the  Crout  LU  factorization  and  use  back  substitution  to  solve  the  system  at  each  time  step.  For  tridiagonal 
matrices  the  factorization  and  the  back  substitution  are  both  order  O(N). 

Remark  2  One  point  to  consider  with  this  method  is  that  if  At  is  too  large,  or  too  small,  we  obtain  instability 
in  the  method  which  can  result  in  noisy  oscillations  before,  or  after,  the  actual  signal.  In  our  testing  we  have 
determined  that  there  is  an  optimal  range  of  values  for  At  given  a  fixed  h. 

The  inverse  problem  has  difficulties  if  the  sample  time  of  the  given  data  (which  is  st  :=  1/sr,  where  sr  is 
the  sample  rate  with  units  of  number  of  observations  per  second)  is  smaller  than  the  time  step.  In  this  case 
we  discard  the  extra  data  points,  effectively  multiplying  the  sample  rate  by  an  integer  multiple.  (We  refer  to 
this  new  sample  rate  as  the  effective  sample  rate  sre  which  will  be  used  to  determine  the  observations  of  the 
simulations.)  Note  that  since  dt  =  0(l/iV),  a  larger  N  may  result  in  more  data  points  being  retained,  and 
thus  provides  an  additional  advantage  in  finding  an  optimal  solution. 


3.3  Method  of  Mappings 


We  want  the  discretization  of  the  domain  to  incorporate  the  interfaces  between  the  material  and  the  vacuum 
so  that  physical  characteristics  are  not  split  across  more  than  one  finite  element.  However,  we  do  not  want 
to  limit  ourselves  to  simulating  cracks  that  are  in  effect  multiples  of  a  fixed  h.  Decreasing  h  to  compensate 
would  increase  computational  time  unnecessarily.  Instead  we  employ  the  “method  of  mappings” .  In  essence 
we  map  our  computational  domain  to  a  reference  domain  with  a  predetermined  discretization,  using  a  linear 
transformation.  In  this  way  we  can  ensure  that  all  of  our  boundaries  are  mapped  to  node  locations.  Also, 
regardless  of  how  small  the  crack  may  be,  we  are  guaranteed  a  fixed  number  of  elements  to  accurately  resolve 
the  behavior  of  the  electric  field  in  that  region.  Each  interval  [2j,  Sj-i-i]  is  mapped  to  a  predeterminined  interval 
[zi,zi+ 1]  by 


2  =  zi+ 1 


2  —  Zj 


Zj+1  -  Zj 


+  Zj 


z  -  Zj+ 1 
Zj  -  Zj+1 


In  effect  we  are  using  a  variable  mesh  size,  thus  we  must  take  care  to  check  that  in  mapping  our  domain 
we  do  not  “stretch”  or  “shrink”  any  interval  too  much  because  our  error  is  based  on  the  mesh  size  in  the 
original  domain.  Although  our  reference  domain  may  be  equally  spaced  with  small  subintervals,  in  general 
our  effective  mesh  sizes  can  be  quite  large.  Therefore,  we  monitor  the  magnitude  of  the  scaling  factors  and 
rediscretize  our  reference  domain  if  any  mesh  size  is  too  large.  The  method  of  mappings  also  allows  us  to 
easily  normalize  our  domain  length,  which  justifies  our  theoretical  development  of  the  problem  in  the  domain 
z  e  [0,  !]• 

The  method  of  maps  also  alters  our  inner  products  used  in  the  weak  formuation  of  the  problem.  We 
subsequently  use  a  weighted  inner  product  defined  as  follows: 


C25 

<</>,  V>>=  /  (j)(z)i/j(z)dz 

J  z0 


l 


where  f(z)  =  z  is  the  piece-wise  linear  transformation  on  [20,25]- 
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3.4  Numerical  Simulations 


The  following  figures  depict  the  numerical  solution  of  the  amplitude  of  the  electric  field  at  various  times 
(Figures  6  and  7),  as  well  as  the  signal  recorded  at  receivers  located  at  z  =  0  (Figure  8)  and  z  =  1  (Figure 
9).  We  considered  a  Debye  medium  with  the  following  parameters: 

78.2, 

5.5, 

1  x  1CT5, 

3.16  x  1(T8, 

2  GHz, 

and  a  material  gap  located  at  [23,  24]  =  [.6,  .61]. 

We  have  used  2  GHz  simply  so  that  our  computational  domain  of  z  €  [0, 1]  would  not  have  to  be  scaled 
for  this  demonstration.  Also,  in  practice,  one  would  not  compute  a  domain  so  much  larger  than  the  material, 
just  as  in  an  experiment  the  sensors  should  be  as  close  to  the  material  as  possible  to  reduce  noise.  We  did 
so  here  merely  so  that  the  full  wavelength  of  the  signal  would  be  visible.  Finally,  we  have  carried  the  time 
length  well  past  the  first  reflections,  again  just  for  demonstration  purposes.  We  chose  to  plot  Figures  8  and 
9  with  the  time  unit  scaled  seconds  (ss)  which  is  just  seconds  multiplied  by  the  speed  of  light  in  a  vacuum, 
ss  =  c  ■  s.  This,  along  with  the  placement  of  the  material  .5m  from  the  sensor,  validates  the  timing  of  the 
first  reflection  (see  Figure  8)  in  that  it  has  traveled  lm  in  c  seconds. 


4  Problem  1 


We  now  apply  an  optimization  routine  to  the  least  squares  error  between  a  simulated  signal  and  the  given  data 
to  try  to  determine  the  crack  characteristics.  In  particular  we  will  be  trying  to  find  the  depth,  d  :=  Z2  —  z\, 
and  the  width,  S  :=  z-$  —  Z2,  which  will  produce  a  simulated  signal  most  closely  similar  (in  the  least  squares 
sense)  to  the  data. 

4.1  Inverse  Problem 

All  of  the  following  are  solved  with  respect  to  a  reference  problem  (i?l)  with  these  parameter  values  (see 
also  Figure  4): 


zo  =  0,  z\  =  .2,  Z2  =  -3,  Z3  =  .5, 24  =  .8,  z5  =  1.0 
/  =  4  GHz,  tf  is  one  period 

r  =  8.1  x  10— 12 ,  a  =  1  x  10-5,  es  =  80.1,  =  5.5  in  the  material 

N  =  1024,  Nt  =  12926 

The  sample  rate  for  the  data  is  one  sample  per  .05ns 

With  this  choice  of  parameters,  the  forward  solve  solution  at  z  =  0  clearly  shows  distinct  reflections  from 
the  Zi,  Z2,  and  z 3.  This  clear  distinction  will  aid  in  our  approximating  the  initial  guesses,  thus  making  this 
a  relatively  easy  sample  problem. 

4.1.1  Initial  Guesses 

Assuming  the  physical  parameters  are  given  (either  known  or  from  a  previous  estimation),  we  want  to 
determine  the  depth  of  any  crack  ( d ),  and  the  width  of  that  crack  (6),  using  reflection  and/or  transmission 
signals.  First  we  will  attempt  to  get  very  close  approximations  to  d  and  using  information  about  the  travel 
times  of  the  data  signal,  then  we  will  use  these  values  as  initial  guesses  in  the  optimization  routine. 

We  now  describe  how  to  identify  the  return  time  of  each  reflection,  which  in  turn  is  used  to  estimate 
each  Zj  value,  and  therefore  d  and  d.  We  need  to  be  able  to  identify  the  point  in  time  when  a  reflection 
begins,  i.e. ,  the  root  immediately  in  front  of  a  peak  (or  trough)  in  the  signal.  (Note  that  certain  reflections 
are  oriented  peak  first,  while  others  are  oriented  trough  first.)  Three  methods  of  doing  this  are  as  follows: 

i.  tri(i)  —  mint  such  that  \E(t)\  >  C\  and  t  >  C2,  where  C\  is  chosen  to  avoid  noise  being  confused  with 
reflections,  and  C2  is  large  enough  that  the  previous  reflection  has  passed. 

ii.  trzii)  =  maxf  such  that  \E(t)\  <  C\  and  t  <  C2,  where  C\  is  chosen  to  avoid  noise  being  confused 
with  reflections,  and  C2  is  small  enough  that  the  peak  (or  trough)  has  not  yet  passed. 

iii.  tr^i)  =  max tj  such  that  E(tj)E(tj+ 1)  <  0  and  tj  <  C2,  where  C2  is  small  enough  that  the  peak  (or 
trough)  has  not  yet  passed. 

In  our  implementation  we  chose  a  combination  of  these  methods  in  that  we  take  the  largest  t  that  satisfies 
at  least  one  of  the  above  conditions.  This  seems  to  work  better  than  any  one  method  individually,  especially 
in  the  presence  of  minor  noise.  Note  also  that  certain  reflections  are  oriented  peak  first,  while  others  are 
trough  first,  thus  in  some  cases  C2  is  small  enough  that  the  trough  has  not  yet  passed. 

We  then  use  the  difference  between  the  return  times  of  two  reflections,  multiplied  by  the  speed  of  light 
in  the  medium  (or  vacuum,  appropriately),  to  determine  the  additional  distance  travelled  by  the  latter 
reflection.  Namely,  this  is  twice  the  distance  between  the  two  interfaces.  In  this  manner  we  estimate  d  and, 
when  possible,  5. 

For  small  <5  the  two  reflections  are  overlapping  (see,  for  instance,  Figure  8).  In  these  cases  we  approximated 
the  location  of  the  start  of  the  second  reflection  by  the  location  of  the  trough.  While  this  is  a  rough 
approximation,  for  this  problem  it  was  sufficient.  In  Problem  2  we  will  develop  more  sophisticated  methods 
which  could  also  be  applied  here. 
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4.1.2  Optimization  of  Least  Squares  Error 


With  these  initial  estimates  to  d  and  S  established,  we  now  define  our  inverse  problem  to  be:  find  q  =  { d ,  5}  £ 
Qad  such  that  the  following  least  squares  error  between  the  simulation  and  the  observed  data  is  minimized: 

1  m  s 

=  (6) 

1=1 *= 1 


Here  the  Ejj  are  measurements  of  the  electric  field  taken  at  M  specific  locations  (e.g.  z®  =  0  and/or  z§  =  1) 
and  S  distinct  times  (e.g.,  every  sre  =  0.06ps).  The  E{ti,z°-,q)  are  solutions  of  the  simulations  evaluated 
at  the  same  locations  and  times  corresponding  to  the  data,  El3 ,  and  using  parameter  values  q.  The  set  Qad 
is  the  feasible  set  of  q  values  determined  such  that  d  and  S  are  realistic  (e.g.  positive). 

We  apply  an  inexact  Gauss-Newton  iterative  method  to  the  optimization  problem.  That  is,  we  re-write 
the  objective  function  as 


J(q) 


l 

AES 


rtr 


where  Rk  =  ( E(ti ,  Zf;  q)  —  Eij)  for  k  =  i  +  (j  —  1  )M  is  the  residual.  To  update  our  approximation  to  q  we 
make  the  Inexact  Newton  update  step  q+  =  qc  +  sc  where 


Sc  =  -{R\qc)rR'(qc))  'Vjy 
=  -  (R' (qc)T R' (qc))1  R1  (qc)T  R(qc) 


is  the  step,  qc  is  the  current  approximation,  and  q+  is  the  resulting  approximation.  This  is  an  inexact  method 
because  we  have  disregarded  the  S  Hessians  of  ( E(ti,Zj\q )  —  £)),  which  is  generally  acceptable  for  small 
residual  problems  [7]. 

In  this  simple  case  we  have  a  2  x  2  matrix  inverse,  so  we  can  compute  it  explicitly.  Each  iteration  requires 
one  function  evaluation  and  a  forward  difference  gradient,  which  is  two  additional  function  evaluations  (since 
we  have  two  parameters).  Each  function  evaluation  is  equivalent  to  a  simulation.  Therefore  we  want  as  few 
iterations  as  possible. 


4.1.3  Convergence 

Initial  testing  (with  2  =  0  data  only)  shows  convergence  to  8  decimal  places  of  each  parameter  in  6  iterations 
of  Gauss-Newton  with  initial  guesses  having  at  most  5%  relative  error  in  6  and  2%  in  d.  The  algorithm  does 
not  converge  to  the  correct  solution  if  the  initial  guess  for  d  has  5%  relative  error  or  if  the  initial  guess  for  6 
has  10%  relative  error. 

One  reason  that  the  algorithm  fails  to  converge  is  that  this  objective  function  is  poorly  behaved.  In 
Figures  10  and  11  we  show  plots  of  the  objective  function  with  respect  to  6.  The  two  very  large  peaks  in 
J  on  either  side  of  the  exact  minimizer  are  due  to  the  simulated  solution  going  in  and  out  of  phase  with 
the  exact  solution.  For  this  example,  the  simulated  solution  is  most  out  of  phase  with  the  exact  solution  at 
S  =  .181  and  S  =  .219  which  correspond  to  the  first  and  second  peak  in  J  respectively.  The  signal  received 
at  2  =  0  for  these  two  S  values  are  plotted  versus  the  actual  signal  with  6  =  .2  in  Figures  12  and  13.  The 
same  phenomenon  occurs  in  the  d  direction,  for  the  same  reasons. 

One  may  consider  adding  the  information  at  the  2  =  1  signal  to  the  least  squares  objective  function 
hoping  to  lessen  the  effects  of  the  peaks  in  J.  But  as  demonstrated  in  Figures  14  and  15,  which  show  the 
objective  function  using  both  signals  and  using  just  the  2  =  1  signal  respectively,  the  out  of  phase  phenomena 
is  still  present. 

Very  few  optimization  routines  can  expect  global  convergence  without  having  initial  conditions  sufficiently 
in  between  the  two  peaks  in  J.  The  effective  convergence  region  for  this  objective  function  applied  to  this 
problem  (with  or  without  observations  at  2  =  1)  is  within  about  8%  of  the  actual  value  of  S  when  d  is  exact 
and  within  about  7.5%  of  the  actual  value  of  d  when  6  is  exact. 

Note  also  that  the  convergence  region  is  very  dependent  on  the  frequency  of  the  interrogating  signal;  for 
higher  frequencies,  the  region  is  even  smaller.  This  is  because  the  distance  between  the  two  peaks  in  J  is 
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z2e  =  0.3, ij  =  1,  N=1024,  Nt=1 5361 ,  Ns=215,  MinJ(R)  =  0  <§>  delta  =  0.2 


Figure  10:  Nonlinear  Least  Squares  objective  function  versus  5  for  a  large  range  of  S  values  (data  at  z  =  0 
only) 


z2e  =  0.3, iJ  =  1,  N=1024,  Nt=1 5361 ,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  11:  Nonlinear  Least  Squares  objective  function  versus  S  for  a  small  range  of  6  values  (data  at  z  =  0 
only) 


linearly  dependent  on  A,  the  wave  length  of  the  interrogating  signal.  This  has  profound  implications  for  our 
desire  to  interrogate  with  signals  in  the  terahertz  range. 

Further,  this  is  considering  only  a  one  parameter  minimization  problem  with  the  other  parameter  held 
fixed  at  the  exact  solution.  Convergence  is  much  worse  for  the  actual  two  parameter  problem.  The  difficulties 
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Comparison  of  si  for  different  deltas  (using  215  sample  points) 


Figure  12:  Signal  received  at  z  =  0  for  two  different  values  of  S  demonstrating  the  simulation  going  out  of 
phase  on  the  left  with  the  data  signal 


Comparison  of  si  for  different  deltas  (using  215  sample  points) 


Figure  13:  Signal  received  at  z  =  0  for  two  different  values  of  6  demonstrating  the  simulation  going  out  of 
phase  on  the  right  with  the  data  signal 


with  the  least  squares  objective  function  in  this  problem  are  clearly  visible  in  the  surface  plots  shown  in 
Figures  16  and  17.  This  explains  why  d  has  a  smaller  convergence  region,  since  J  is  much  more  sensitive  to 
d.  Also  we  see  that  there  are  peaks  in  J  in  the  d  direction  as  well. 
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z2e  =  0.3, ij  =  1 ,02  =  1 ,  N=1 024,  Nt=1 5361 ,  Ns=21 5,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  14:  Nonlinear  Least  Squares  objective  function,  using  signals  received  at  both  z  =  0  and  z  =  1,  versus 
6 


z2e  =  0.3, iJ  =  1,01  =  0,02  =  1,  N=1024,  Nt=15361,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  15:  Nonlinear  Least  Squares  objective  function,  using  signals  received  at  z  =  1  only,  versus  S 
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delta 


Figure  16:  Nonlinear  Least  Squares  objective  function  (J)  versus  S  and  d  projected  to  2D  (data  at  z  =  0 
only) 


Figure  17:  Nonlinear  Least  Squares  objective  function  (J)  versus  <5  and  d  in  3 D  (data  at  z  =  0  only) 
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The  diagonal  “trench”  occurs  approximately  along  the  line 

d=-l-(6-6*)+d*  (7) 

(where  *  denotes  the  exact  solution).  This  happens  because  when  the  depth  d  is  decreased  by,  say  a,  the  2 1 
and  22  reflections  of  the  simulation  no  longer  correspond  with  those  of  the  data  (think  of  the  two  overlapping 
reflections  at  t  =  1.5ss  in  Figure  8  being  shifted  slightly  to  the  left  in  some  simulation). 

With  d  decreased  by  a  we  must  increase  S  to  try  to  compensate  and  get  as  much  of  the  signal  from  the 
simulation  to  correspond  with  the  data.  Since  6  only  affects  the  22  reflection,  it  is  clear  that  the  21  reflections 
will  not  correspond.  If  <5  is  increased  so  that  the  distance  traveled  by  the  22  reflection  is  increased  by  6, 
then  we  will  need  for  the  time  required  for  a  signal  to  travel  a  distance  a  in  the  material  to  be  equal  to  the 
time  required  for  a  signal  to  travel  a  distance  b  in  a  vacuum.  Since  the  distance  traveled  by  the  22  reflection 
is  twice  that  of  the  increase  in  6,  we  compute  the  increase  in  S  to  be  In  our  case  =  5.5,  so  the 

increase  is  2.3452a.  This  explains  the  slope  of  approximately  —  A  in  the  relationship  (7). 

4.2  An  Improved  Objective  Function 

As  demonstrated  above,  the  usual  Nonlinear  Least  Squares  objective  function  has  two  large  peaks  in  J  on 
either  side  of  the  exact  minimizer.  The  reason  for  these  peaks  in  J  is  that  the  simulated  solution  goes  in 
and  out  of  phase  with  the  data  as  d  or  d  change.  When  they  are  precisely  out  of  phase,  there  is  a  very  large 
absolute  error,  which  when  squared,  causes  the  objective  function  to  have  large  peaks.  Therefore,  a  logical 
solution  would  be  not  to  consider  the  absolute  error,  but  the  error  of  the  absolute  values,  i.e. ,  the  following 
objective  function: 

Mq)  =  7^Y,{  |£(ti,*?;?)l-|£il)2-  (8) 

°  2=1 

This  will  prevent  the  fact  that  the  signals  go  in  and  out  of  phase  with  each  other  from  having  an  impact  on 
the  objective  function,  since  positive  magnitudes  cannot  cancel  each  other  out.  This  gives  a  more  accurate 
measure  of  the  difference  between  two  signals.  Note  that  A (9)  is  not  differentiable  on  a  set  of  measure  zero; 
this  is  very  unlikely  to  affect  the  finite  difference  computations  of  the  gradients,  and  did  not  seem  to  present 
problems  in  our  testing.  We  plot  J'liq)  versus  S  in  Figures  18  and  19. 

Note  that  while  we  have  effectively  eliminated  the  peaks  on  either  side  of  the  exact  solutions,  in  essence 
we  have  merely  converted  them  to  local  minima!  But,  since  the  minima  occur  for  the  same  reasons  the 
peaks  in  J  had  been  occurring,  they  occur  at  the  same  values  of  S.  Note  that  we  can  see  from  plotting  the 
signals  that  they  were  exactly  out  of  phase  when  they  were  shifted  by  j,  where  A  is  the  wavelength  of  the 
interrogating  signal.  Therefore  <5  is  off  by  4 .  Since  we  determine  the  frequency  of  the  interrogating  signal, 
this  is  a  known  quantity,  and  we  can  predict  where  these  local  minima  will  occur  a  priori! 

Most  optimization  routines  will  continue  until  they  find  a  local  minimum,  and  since  the  two  false  minima 
described  above  are  at  least  close  to  “predictable”  locations,  we  can  easily  test  on  either  side  of  any  detected 
minimum  to  determine  if  it  is  in  fact  global.  If  J2  is  less  at  a  fixed  distance  in  the  6  direction  (e.g.,  |j  on 
either  side  of  a  detected  minimum,  i.e.,  at  either  test  location  of  the  local  minimum,  we  restart  Gauss-Newton 
at  the  new  best  guess.  Thus  if  either  of  the  two  local  minima  described  above  is  found,  we  will  eventually 
have  global  convergence  if  one  of  their  test  locations  are  sufficiently  close  to  the  global  minimizer. 

To  graphically  demonstrate  this  approach,  we  have  added  dotted  and  dash-dotted  lines  to  the  the  graph 
of  .J2(S)  in  Figure  19  to  represent  the  test  values  of  the  first  and  second  local  minima  respectively.  One  can 
clearly  see  that  if  either  local  minimum  is  detected,  the  test  value  either  to  its  right  if  it  is  the  first  one,  or 
to  its  left  if  it  is  the  second,  will  give  a  smaller  J2  and  should  eventually  lead  to  the  global  minimizer  being 
detected.  Using  this  method  we  have  in  principle  increased  our  convergence  region  to  about  25%  of  S  when 
d  is  exact. 

Our  objective  function  also  has  local  minima  in  the  d  direction,  as  demonstrated  in  Figure  20.  Again 
they  are  in  predictable  locations,  being  at  roughly  d*  ±  — ,  where  d*  represents  the  global  minimizer. 

Thus,  if  we  know  A  and  Coo  (which  we  should),  then  we  can  determine  exactly  where  to  test  for  the  global 
minima  once  a  local  minimum  has  been  found.  Using  the  above  method,  we  have  in  principle  increased  the 
d  convergence  region  from  about  7.5%  to  about  15%  when  6  is  exact. 
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d  =  0.1 ,  N=1 024,  Nt=1 5361 ,  Ns=21 5,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  18:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)  versus  S  for  a  large  range  of  S 
values 


d  =  0.1,  N=1024,  Nt=  15361 ,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  19:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)  versus  S  for  a  small  range  of  6 
values.  The  dotted  lines  represent  the  delta  values  that  will  be  tested  if  a  local  minimum  is  found 


Other  possible  modifications  to  the  Least  Squares  objective  function  having  similiar  effects  include  squar¬ 
ing  the  signal  instead  of  taking  the  absolute  value  (thus  preserving  smoothness  everywhere),  or  just  halving 
the  reference  signal  so  that  we  only  have  a  positive  amplitude  to  begin  with. 

In  addition  to  the  two  new  minima  created  above,  we  still  have  other  local  minima  to  deal  with  which 
are  present  in  both  formulations  of  the  objective  function.  Some  of  these  local  minima  can  be  reduced  by 
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N=1024,  Nt= 15361 ,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  20:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)  versus  d  for  a  small  range  of  d 
values.  The  dotted  lines  represent  the  delta  values  that  will  be  tested  if  a  local  minimum  is  found 


using  a  smaller  mesh  size  (more  finite  elements)  since  they  are  the  result  of  numerical  error  (see  Figures  21 
and  22). 

All  of  the  above  discussions  and  Figures  18  through  22  involve  J2  with  observations  only  at  z  =  0.  We 
have  also  done  tests  using  only  the  z  =  1  transmitted  data  in  the  objective  function  from  (8)  (see  Figure 
23) .  Clearly  this  function  does  not  have  the  same  problem  of  the  third  reflections  (since  z  =  1  only  sees  the 
even  numbered  reflections).  But  it  does  still  have  a  very  deep  and  significantly  close  local  minima  of  its  own. 
This  minimum  is  caused  by  the  same  type  of  problem  that  causes  the  z  =  0  minimum.  However,  we  have 
found  that  by  combining  both  z  =  0  and  z  =  1  data  in  the  objective  function  J2  we  can  get  much  better 
results  since  their  deficiencies  somewhat  average  out  (see  Figure  24). 
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z2e  =  0.3,  N=2048,  Nt=30721 ,  Ns=21 5,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  21:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2),  using  a  simulation  with  twice  as 
many  meshes,  versus  S  for  a  large  range  of  S  values 


z2e  =  0.3,  N=2048,  Nt=30721 ,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


delta 


Figure  22:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2),  using  a  simulation  with  twice  as 
many  meshes,  versus  6  for  a  small  range  of  5  values.  The  dotted  lines  represent  the  delta  values  that  will  be 
tested  if  a  local  minimum  is  found 
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z2e  =  0.3,01  =  0,02  =  1 ,  N=1 024,  Nt=1 5361 ,  Ns=21 5,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  23:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)>  using  the  signal  received  at  z  =  1 
only,  versus  S  for  a  large  range  of  6  values 


02=1,  N=1 024,  Nt=1 5361 ,  Ns=21 51 ,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  24:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2),  using  the  signal  received  at  z  =  0 
and  2  =  1,  versus  S  for  a  large  range  of  S  values 
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4.3  Testing  J2 

In  order  to  determine  the  limitations  of  an  optimization  routine  to  minimize  our  objective  function  J2  in 
a  more  practical  setting  we  examine  J2  versus  q  when  error  is  present.  In  particular  we  try  both  adding 
random  noise  to  the  data  signal,  as  well  as  testing  bad  initial  guesses  for  6  and  d.  It  should  be  noted  that 
although  we  have  previously  established  that  there  are  certain  benefits  to  having  data  at  z  =  1  in  the  tests 
reported  on  below  we  assume  that  it  is  not  available  and  used  only  observations  at  z  =  0. 

4.3.1  Sensitivity  to  Initial  Guesses 

For  J2  described  in  (8)  we  compute  a  surface  plot  for  various  6  and  d  values,  which  is  shown  in  Figures  25 
and  26.  We  notice  immediately  that  the  objective  function  is  much  more  sensitive  to  d  than  to  S,  therefore 
it  is  imperative  that  our  initial  guess  for  d  is  as  good  as  possible. 

To  give  an  idea  of  what  may  happen  if  our  d  guess  were  not  within  the  15%  our  testing  has  determined 
is  necessary,  we  plot  the  objective  function  versus  <5  for  three  values  of  d,  which  are  3%,  15%,  and  30%  off 
respectively,  in  Figures  27,  28,  and  29. 

It  is  clear  that  the  15%  case  should  have  no  trouble  converging  with  relatively  good  initial  guess  for  <5, 
but  even  if  the  initial  value  for  S  were  exact,  the  30%  case  could  quite  possibly  converge  to  the  minimum  at 
the  far  left  of  Figure  29.  It  may  be  surprising  that  there  should  be  a  minimum  for  very  small  S  values  at  all, 
even  more  so  that  it  is  in  fact  the  global  minimum!  This  can  be  explained  by  first  remembering  that,  for 
example  in  the  15%  case,  the  original  data  was  computed  with  d  —  .1  but  the  simulations  used  d  =  .085,  so 
we  should  expect  that  regardless  of  <5  there  is  a  part  of  the  simulated  signal  that  will  not  match  the  data, 
namely  the  first  reflection.  Thus  the  first  reflection  of  the  data  is  not  matched  by  the  simulation,  unless  S 
is  small  enough  to  match  it,  e.g.  5  =  .035,  and  this  is  exactly  what  is  happening,  as  displayed  in  Figure 
30.  Figure  31  shows  the  other  local  minimum  of  the  15%  case  where  the  z 3  reflection  of  the  simulation  does 
match  the  Z3  reflection  of  the  data. 

Notice  that  the  distance  between  the  two  minimum  values  of  S  (see  for  example  Figure  28  or  29)  is  exactly 
S*  =  0.2,  which  is  what  would  be  expected.  However,  we  cannot  apply  the  same  idea  as  before  where  we 
add  or  subtract  a  fixed  amount  to  test  for  other  local  minima,  since  for  one,  the  “more  optimal”  of  the  two 
is  farther  from  the  “true”  solution,  and  also,  we  would  have  to  know  S  in  the  first  place  to  add  or  subtract 
it,  but  that  is  what  we  are  trying  to  estimate! 
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Figure  26:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)  versus  S  and  d  in  3 D 
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Comparison  of  si  for  different  deltas  (using  215  sample  points) 


Figure  30:  Signals  received  at  z  =  0  (where  (d,S)  =  (.1,  .2)  corresponds  to  the  given  data  and  (d,  S)  = 
(.085,  .035)  corresponds  to  the  simulation)  demonstrating  the  z<i  reflection  of  the  simulation  corresponding 
to  the  zi  reflection  of  the  data 


Comparison  of  si  for  different  deltas  (using  215  sample  points) 


Figure  31:  Signals  received  at  z  =  0  (where  (d,S)  =  (.1,  .2)  corresponds  to  the  given  data  and  (d,6)  = 
(.085,  .235)  corresponds  to  the  simulation)  demonstrating  the  Zi  reflection  of  the  simulation  matching  the  22 
reflection  of  the  data  and  ignoring  the  z\  reflection 
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4.3.2  Random  Observation  Noise 


In  order  to  test  the  feasibility  of  this  procedure  as  an  estimation  method,  we  have  produced  synthetic  data 
for  our  observations  £%  In  an  actual  experiment,  one  must  assume  that  the  measurements  are  not  exact.  To 
simulate  this  we  have  added  random  noise  to  the  original  signal.  The  absolute  value  of  the  noise  is  relative 
to  the  size  of  the  signal.  If  E,  is  the  data  sampled,  then  we  define  E.t  =  £j(l  +  vrji) ,  where  rii  are  independent 
normally  distributed  random  variables  with  mean  zero  and  variance  one.  The  coefficient  v  determines  the 
relative  magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  in  particular,  v  =  0.05  corresponds 
to  10%  noise  and  v  =  0.025  to  5%  noise. 

Plots  of  the  resulting  objective  functions  for  various  values  of  v  ranging  from  2%  to  40%  are  shown  in 
Figures  32-37.  We  notice  that  the  structure  of  the  curves  is  not  significantly  affected,  nor  is  the  location 
of  the  global  minimum.  However  the  magnitude  of  the  objective  function  is  increased,  making  Inexact 
Newton  methods  slightly  less  reliable  due  to  the  larger  residual.  Still,  our  results  show  that  the  minima  were 
consistently  found  and  within  a  reasonable  amount  of  time.  Select  examples  are  summarized  in  Table  1. 


V 

d 

S 

J 

Iterations 

CPU  Time  (s) 

0 

0.1 

0.2 

1.32319E-10 

7 

160 

0.01 

0.099994 

0.199969 

0.00792792 

8 

186 

0.05 

0.099974 

0.199835 

0.199489 

13 

291 

0.2 

0.099928 

0.199204 

3.04619 

20 

435 

Table  1:  Number  of  Iterations  and  CPU  Time  for  Gauss-Newton  given  various  relative  magnitudes  of  random 
error 
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z2e  =  0.3, nu  =  0.01,  N=1024,  Nt=15361,  Ns=215,  MinJ{R)  =  0.0280271  @  delta  =  0.2 


Figure  32:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2),  with  2%  random  noise,  versus 
for  a  large  range  of  S  values 


z2e  =  0.3, nu  =  0.01 ,  N=1 024,  Nt=1 5361 ,  Ns=21 5,  MinJ(R)  =  0.0280271  @  delta  =  0.2 


delta 

Figure  33:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2),  with  2%  random  noise,  versus  5 
for  a  small  range  of  6  values 
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Figure  38:  Schematic  of  Problem  2:  determining  the  depth  and  width  of  a  gap  between  a  dielectric  slab  and 
a  metallic  backing 


5  Problem  2 

We  will  apply  the  most  useful  techniques  from  Problem  1  to  a  new  formulation  of  the  interrogation  problem. 
In  Problem  2  we  consider  a  dielectric  slab  and  a  metallic  backing  (conductor)  with  a  possible  gap  between 
the  two  (see  Figures  38  and  39).  Applications  of  this  specific  formulation  included  detecting  delamination 
of  insulation  from  metallic  containers,  e.g.,  insulating  foam  on  a  space  shuttle  fuel  tank.  In  order  for  this 
numerical  approach  to  be  useful  in  this  particular  application  we  must  be  able  to  resolve  a  slab  thickness  of 
at  least  20cm  and  a  frequency  of  100 GHz. 

We  will  again  assume  the  same  physical  parameters  for  our  dielectric  and  consider  the  gap  as  a  vacuum. 
The  variables  d  and  <5  are  still  the  depth  and  the  width  of  the  gap  respectively.  One  major  difference  is  that 
in  this  problem  we  are  only  able  to  detect  the  electromagnetic  signal  in  front  of  the  material.  Also,  since  the 
metallic  backing  reflects  much  of  the  signal,  we  have  considerably  more  overlapping  of  the  reflections  to  worry 
about.  These  properties  contribute  to  the  fact  that  this  formulation  leads  to  a  much  more  difficult  inverse 
problem.  For  this  reason  we  will  be  using  more  sophisticated  optimization  routines  including  a  Levenberg- 
Marquardt  parameter  and  Implicit  Filtering.  We  will  also  need  to  develop  different  approximation  methods 
for  our  initial  guesses. 

The  implementation  of  this  problem  has  several  minor  differences  from  the  previous  one.  First,  we 
now  only  need  to  represent  two  interfaces  Z\  and  z 2 ,  with  zg  and  £3  being  the  front  and  back  boundaries, 
respectively.  Thus  now  we  define  the  depth  of  the  gap  as  d  :=  z^  —  Z\  and  the  width  as  (5  :=  £3  —  £2-  Also, 
as  previously  mentioned,  the  conductive  metal  backing  reflects  the  signal,  and  hence  we  must  change  our 
absorbing  boundary  conditions  at  z  =  1  (for  a  finite  computational  domain),  to  an  actual  fixed,  Dirichlet 
boundary  condition  (E  =  0).  We  must  modify  our  finite  element  matrices  accordingly,  as  well.  Otherwise, 
the  numerical  method  for  simulation  is  the  same  as  it  was  for  Problem  1,  namely  standard  finite  element 
methods  for  spatial  derivatives,  and  an  alternating  implicit /explicit  centered  difference  time  stepping  scheme. 
Sample  solutions  are  plotted  in  Figure  40. 

We  again  define  our  inverse  problem  to  be:  find  q  :=  {d,  5}  £  Qad  such  that  an  objective  function 
representing  the  error  between  the  simulation  and  the  observed  data  is  minimized: 

min  J(q). 

ad 

Here  the  measurements  of  the  electric  held,  £),  are  taken  only  at  z  =  0,  but  still  at  S  distinct  times  (e.g., 
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0  0.005  0.01  0.015  0.02 

z  (meters) 


Figure  39:  The  domain  of  the  material  slab  with  a  gap  between  the  medium  and  a  metallic  conductive 
backing:  12  =  {z\zi  <  z  <  z-i}  ■ 


every  0.06ps).  The  solutions  of  the  simulations,  E(ti,0\q),  are  evaluated  at  the  same  location  and  times 
corresponding  to  the  given  data,  and  using  parameter  values  q.  In  lieu  of  actual  data  from  experiments,  we 
again  create  our  observed  data  by  using  the  simulator,  however,  the  only  information  that  is  given  to  the 
minimizer  is  the  data  observed  at  z  =  0,  which  we  will  denote  by  E. 

The  system  that  we  use  to  represent  the  propagation  of  the  electric  field,  and  thus  simulate  in  order  to 
solve  our  inverse  problem,  is  as  follows,  and  includes  the  above  mentioned  Dirichlet  condition  at  z  =  1: 

/x0eo(l  +  (eoo  -  1  )Iq)E  +  no IqP  +  Moc-E  -  E"  =  -p,0  js 

tP  +  P  =  re0(es  -  e^E 
[E  -  cE']z-o  =  0 
[E\z= l  =  0 

E(0,z)  =  0 
E(0,z)  =  0 

with 

Js(t,z)  =  6(z)sin(ojt)I[0,tf](t). 

See  Section  2  for  a  complete  description. 

5.1  Objective  Function 

As  in  the  previous  problem,  we  encounter  difficulties  when  attempting  to  use  the  standard  Least  Squares 
objective  function  to  compute  the  error  between  the  simulated  signal  and  the  observed  data.  The  constructive 
interference  of  peaks  and  troughs  produces  peaks  in  J  in  the  objective  function  on  all  sides  of  the  global 
minimum  which  make  it  nearly  impossible  to  find  the  solution  in  the  middle.  The  peaks  in  J  in  the  d 
direction  are  clearly  apparent  in  Figure  41.  In  contrast,  Figure  42  shows  a  surface  plot  of  our  modified  least 
squares  objective  function 

g 

Mq)  =  7^^2\\E(ti,0;q)\-\Ei\  . 

i=  1 

Remark  3  The  gradual  slope  down  and  to  the  left,  in  the  surface  plots  is  due  to  the  increased  attenuation  of 
signals  as  the  depth  is  increased.  Also,  the  drop  off  at  large  depth  values  is  an  artifact  of  having  a  finite  end 
time  for  collecting  data,  i.e.,  the  reflection  from  the  back  of  the  material  did  not  return  in  time. 
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Figure  40:  Computed  solutions  at  different  times  of  a  windowed  electromagnetic  pulse  incident  on  a  Debye 
medium  with  a  crack  between  the  medium  and  a  metallic  conductive  backing.  The  width  of  the  slab  is 
d  =  .02m  and  the  width  of  the  gap  is  <5  =  .0002  (barely  visible  at  the  far  right  of  the  gray  region). 


It  is  clear,  as  before,  that  the  initial  guess  is  crucial  to  the  success  of  any  optimization  routine.  To  further 
demonstrate  this  importance,  we  expand  our  surface  plot  out  in  the  5  direction,  and  zoom  in  with  respect 
to  d.  Figures  43  and  44  show  Ji  and  J2  respectively.  Notice  that  although  J2  does  not  exhibit  the  familiar 
peaks  in  J  of  Ji,  it  does  however  still  have  many  local  minima,  which  are  as  difficult,  if  not  more  so,  to  avoid 
in  a  minimization  routine. 

The  local  minima  occur  approximately  every  ^  along  the  line 


d  = 


S  +  b. 


This  happens  for  the  same  reason  as  in  Problem  1  (see  Section  4.1.3).  For  an  illustration,  see  Figures  45  and 
46.  In  each  figure  we  see  that  the  small  objective  function  value  is  due  primarily  to  the  coinciding  of  the 
single  largest  peak  and/or  the  largest  trough.  (Note  that  the  trailing  oscillations  are  not  numerical  noise; 
the  reflections  appear  as  they  do  because  S  «  |  «  7.5  x  10-4.) 

Because  we  cannot  eliminate  these  local  minima,  we  must  appeal  to  the  procedure  that  worked  in  the 
previous  problem,  namely  testing  “check  points”.  Since  we  know  where  these  local  minima  are  occurring 
with  respect  to  the  global  minimum,  if  our  minimization  routine  finds  what  it  suspects  to  be  a  local  minima, 
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J1  depth 


depth 


Figure  41:  Surface  plot  of  Least  Squares  objective  function  demonstrating  peaks  in  J. 
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Figure  42:  Surface  plot  of  modified  Least  Squares  objective  function  demonstrating  lack  of  peaks  in  J. 


say  (di,<5i),  we  simply  check  (di  ±  =F  where  a  =  l/i/l  +  Coo-  If  we  find  a  lower  objective 

function  value,  we  restart  our  optimization  routine  at  that  “check  point”. 
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Figure  43:  Close  up  surface  plot  of  Least  Squares  objective  function  demonstrating  peaks  in  J,  and  exhibiting 
many  local  minima. 


J2 


Figure  44:  Close  up  surface  plot  of  modified  Least  Squares  objective  function  demonstrating  lack  of  peaks 
in  J,  but  exhibiting  many  local  minima. 


5.2  Initial  Guesses 

In  spite  of  our  faith  in  the  “check  point”  method,  we  still  desire  to  find  the  best  initial  guesses  for  our 
optimization  routine  as  possible  so  that  we  may  hopefully  find  the  global  minimum  without  restarting.  As 
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Signal  received  at  z=0 


t  (ns) 

Figure  45:  Data  from  (d*,S*)  and  a  simulation  from  the  “check  point”  (d*  —  a^,6*  +  The  first 

trough  cannot  be  matched,  but  S  is  sufficiently  large  so  that  the  signal’s  peak  matches  with  that  of  the  data. 


before,  we  use  the  travel  time  of  the  first  trough  to  approximate  the  location  of  the  first  interface.  However, 
in  this  formulation  we  can  take  advantage  of  some  of  the  characteristics  of  the  signals.  For  example,  the 
first  reflection  off  of  the  gap  is  always  trough-first,  and  the  second  (as  well  as  each  subsequent  reflection) 
is  always  peak-first.  For  this  reason,  if  we  want  to  locate  the  first  trough  we  can  simply  find  the  largest 
peak  (belonging  to  the  second  reflection)  and  back  track.  It  is  a  very  simple  matter  to  find  a  maximum  or 
minimum  of  a  vector  of  values.  After  the  location  of  the  largest  peak  is  found,  we  back  track  to  find  the 
minimum  in  front  of  it,  namely  that  belonging  to  the  first  reflection  off  of  the  gap.  Then  using  the  procedure 
described  in  Section  4.1,  we  approximate  the  root  immediately  in  front  of  this  trough.  That  gives  us  the 
travel  time  for  the  first  reflection  off  of  the  gap,  which  in  turn  gives  us  the  depth  d  of  the  gap. 

Finding  6  is,  unfortunately,  not  nearly  as  straightforward.  There  are  two  main  possibilities,  and  therefore, 
two  differing  approaches  to  approximating  5: 

(i)  The  leading  trough  of  the  first  reflection  and  the  second  reflection  are  disjoint  (i.e.,  5  >  -|).  In  this 
case  we  can  find  the  locations  of  the  peak  and  trough  and  use  the  travel  time  between  the  two  to 
approximate  6.  We  denote  this  approximation  by  <5i .  (Note  that  the  observed  peak  is  not  necessarily 
the  same  as  the  original  peak  unless  5  >  j,  but  it  is  still  a  good  approximation).  See  Figure  47. 

(ii)  The  second  reflection  partially  truncates  the  trough  of  the  first  (i.e.,  S  <  |).  Asa  rough  approximation, 
we  can  assume  that  the  location  of  the  actual  minimum  (trough)  is  where  the  two  signals  begin  to 
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Figure  46:  Data  from  (d*,  5*)  and  a  simulation  from  the  “check  point”  (d*  —  2a  8*  +  2 Again,  the 
first  trough  cannot  be  matched,  but  this  time  simulated  signal  has  no  cancellations  so  that  its  largest  peak 
matches  with  that  of  the  data. 


interfere  with  each  other  (the  observable  minimum).  See  Figure  48.  We  denote  this  approximation  by 

64. 

A  more  accurate  method  is  to  use  triangles  to  approximate  the  two  reflections.  By  knowing  the 
location  of  the  maximum  and  minimum  (peak  and  trough,  respectively),  and  also  the  beginning  of 
the  first  signal  (from  Section  4.1)  and  the  rough  approximation  to  the  beginning  of  the  second  signal 
using  Si i,  we  can  estimate  the  slopes  of  the  two  triangles  with  finite  differences.  Also  note  that  since 
the  two  signals  are  added,  the  observed  root  between  the  peak  and  trough  in  the  combined  signal  is 
actually  an  equilibrium  point  between  the  two  signals.  By  setting  equal  to  each  other  the  two  linear 
approximations  for  each  of  the  two  signals,  evaluated  at  the  equilibrium  point,  we  can  solve  for  the 
distance  between  the  starting  point  of  each  signal,  and  thus  for  S3.  See  Figure  49.  Specifically,  let 
(pi,qi)  be  the  location  of  the  trough  of  the  combined  signal  and  (p2,<72)  be  the  location  of  the  peak. 
Let  n  be  the  location  of  the  root  in  front  of  the  trough,  and  r2  be  the  root  between  the  trough  and 
peak.  Estimate  the  slope  of  the  first  signal,  mi  <  0,  using  (pi,  q\ )  and  n.  Now  if  we  let  y  =  r2  —  rq 
and  say  x  is  the  actual  distance  between  n  and  the  beginning  of  the  second  signal,  then  setting  the 
linear  approximations  equal  in  magnitude,  but  opposite  in  sign,  at  r2  yields 

-miy  =  m2(y  -  x). 
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Figure  47:  The  top  plot  represents  several  signals  which  may  be  observed  in  a  simulation  of  Problem  2.  The 
bottom  plot  shows  the  sum  of  the  top  signals.  The  peak  of  the  second  signal  is  just  beginning  to  be  obscured 
by  the  first  when  S  becomes  less  than  A/4.  Thus  the  observable  maximum  is  still  a  good  approximation  of 
the  peak  of  the  second  signal,  and  a  trough  to  peak  distance  can  be  used  to  estimate  <5. 


Now  we  can  estimate  the  slope  of  the  second  signal,  m2  >  0,  using  (p2,(72)  and  ( V2 ,  —miy).  Also,  we 
can  re-write  the  above  equation  as 

f  —mi  +  m2\ 

x  =  -  y- 

\  m2  J 

To  find  S3  we  simply  divide  a;  by  2  and  the  (scaled)  speed  of  light  in  the  material,  i.e.,  y/e^. 

Since  each  of  the  two  situations  above  is  dependent  on  the  parameter  it  is  approximating,  we  must  also 
determine  which  of  the  above  methods  is  most  appropriate  to  use.  Thus  we  use  the  most  precise  of  the 
available  methods  to  determine  the  situation,  i.e.,  64,  instead  of  £3  since  in  general  £3  underestimates  S  so  we 
do  not  want  to  use  it  as  a  criterion  for  determining  whether  <5  is  small.  (Note  that  when  <5  is  indeed  small, 
S 3  is  more  accurate  than  ^4.)  The  estimate  for  64  tends  to  be  an  overestimate,  and  is  only  valid  if  S  <  ^ . 
Unfortunately,  also  tends  to  be  an  overestimate,  so  we  prefer  to  only  trust  it  entirely  if  it  is  larger  than 
)) .  If  neither  nor  (>3  is  a  sufficient  approximation  we  choose  to  use  the  average  of  the  two,  and  call  it  <J2. 

Therefore  our  algorithm  for  approximating  <5  is  as  follows: 

(a)  If  64  <  |  then  use  <53 

(b)  else  if  ch  >  ^  then  use  <5i 
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t  (ns) 

Figure  48:  The  top  plot  represents  several  signals  which  may  be  observed  in  a  simulation  of  Problem  2.  The 
bottom  plot  shows  the  sum  of  the  top  signals.  The  trough  of  the  first  signal  is  partially  truncated  by  the 
second  signal.  In  this  case  the  observed  minimum  is  a  still  a  good  approximation  to  where  the  second  signal 
begins.  For  smaller  8,  a  linear  approximation  must  be  used. 


(c)  else  use  82  (average  between  cq  and  S3). 

We  tested  our  approximating  methods  on  exact  depth  (d)  values  of:  .02,  .04,  .08,  .1,  and  .2  to,  and 
values  of  width  (5):  .0001,  .0002,  .0004,  .0006,  and  .0008  to.  Since  |  is  the  transition  point  between  the 
two  situations,  it  is  understandable  why  8  close  to  this  value  is  the  most  difficult  to  accurately  resolve.  We 
chose  this  range  of  S' s  because  our  choice  of  frequency  gives  |  =  3.7475  x  10_4to.  See  Tables  2  and  3  for 
the  initial  estimates  of  d  and  8  respectively. 

The  approximations  improve  slightly  as  the  number  of  finite  elements  is  increased,  and  seemed  to  converge 
to  fixed  values.  This  suggests  that  numerical  error  (and  instability)  can  affect  the  estimates.  For  each  case 
there  is  a  significant  amount  of  numerical  error  in  the  simulations  below  a  certain  number  of  elements, 
therefore  in  approximating  8  we  chose  to  use  the  number  of  elements  just  above  the  threshold.  The  number 
of  finite  elements  ( TV )  used  in  the  simulations  is  given  for  each  case.  Since  approximating  d  does  not  depend 
on  5 ,  we  show  the  progression  of  estimates  as  N  increases.  The  values  in  Table  2  which  are  bold  denote  the 
values  used  to  approximate  5  in  Table  3. 

These  initial  estimates  may  seem  relatively  inaccurate,  some  8  approximations  being  almost  100%  off 
from  the  actual  value,  but  in  these  tests  all  were  sufficiently  close  to  the  global  minima  to  not  result  in  a 
false,  local  minima.  While  our  “check  point”  method  is  available  if  needed,  it  is  much  more  efficient  to  have 
an  accurate  initial  estimate  than  to  restart  after  optimizing  from  a  bad  one.  Figures  50  through  54  display 
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Figure  49:  This  schematic  shows  the  roots,  extrema,  distances,  and  slopes  used  in  the  computation  of  83. 


N 


d 

1024 

2048 

4096 

8192 

16384 

.02 

0.0200036 

0.019996 

0.0199998 

0.0199998 

.04 

0.0399689 

0.0399766 

0.0399919 

0.0399958 

.08 

.1 

.2 

0.0799225 

0.0799609 

0.0799839 

0.0799992 

0.0999799 

0.200006 

Table  2:  The  initial  estimates  of  d.  The  values  in  bold  denote  the  values  used  to  approximate  5  in  Table  3. 

<5 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.41856e-05 

0.000127645 

0.000611592 

0.000683544 

0.000791472 

.04 

(N=2048) 

0.000113806 

0.000210178 

0.000611592 

0.000647568 

0.000791472 

.08 

(N=4096) 

2.89938e-05 

0.00019437 

0.000575616 

0.000647568 

0.000755496 

.1 

(N=8192) 

0.000151182 

0.000239409 

0.000593604 

0.000647568 

0.000791472 

.2 

(N=16384) 

8.74169e-05 

0.000146323 

0.000575616 

0.000647568 

0.000791472 

Table  3:  The  initial  estimates  of  5. 


surface  plots  of  the  modified  Least  Squares  objective  function  zoomed  in  around  the  region  that  our  initial 
estimates  are  located.  It  is  evident  that  relatively  small  errors  in  the  initial  estimates  could  prevent  an 
optimization  routine  from  finding  the  global  minima.  For  example,  in  Figure  50  if  the  initial  estimate  for  d 
were  1  x  10~4  too  high  and  <5  too  high  by  2  x  10-4  then  the  gradient  would  cause  both  8  and  d  to  increase, 
thus  leading  to  the  wrong  minima.  In  cases  such  as  these,  the  “check  point”  method  provides  a  last  resort. 
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d=.02,8=.0001 


Figure  50:  Very  close  surface  plot  of  the  modified  Least  Squares  objective  function  using  d  =  .02  and 
6  =  .0001. 


d=.02,8=.0002 
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Figure  51:  Very  close  surface  plot  of  the  modified  Least  Squares  objective  function  using  d  =  .02  and 
S  =  .0002. 


42 


d=.02,8=.0004 
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Figure  52:  Very  close  surface  plot  of  the  modified  Least  Squares  objective  function  using  d  =  .02  and 
5  =  .0004. 


d=. 02,5=. 0006 


delta 
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Figure  53:  Very  close  surface  plot  of  the  modified  Least  Squares  objective  function  using  d  =  .02  and 
S  =  .0006.  (Note  that  the  axis  is  shifted  from  the  previous  Figures  in  order  to  include  the  minimum.) 
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d=. 02,5=. 0008 


0.0201 


Figure  54:  Very  close  surface  plot  of  the  modified  Least  Squares  objective  function  using  d  =  .02  and 
6  =  .0008.  (Note  that  the  axis  is  shifted  from  the  previous  Figures  in  order  to  include  the  minimum.) 


5.3  Optimization  Method 

Now  that  we  have  approximated  our  initial  guesses,  we  need  to  minimize  the  objective  function  in  order  to 
solve  the  inverse  problem.  In  Problem  1,  Gauss-Newton  was  sufficient  to  find  the  global  minimum  for  most 
cases.  In  this  formulation,  however,  we  will  apply  more  sophisticated  methods,  reverting  to  Gauss-Newton 
whenever  possible  since  its  convergence  rate  is  best. 

The  first  modification  we  make  to  Gauss-Newton  is  to  add  a  Levenburg-Marquardt  parameter,  vc  (see 
[7]).  The  Inexact  Newton  step  becomes 

sc  =  -  (- R'{qc)TR\qc )  +  vc  1)^  R\qc)TR(qc )■ 

The  parameter  adds  regularization  by  making  the  model  Hessian  positive  definite.  The  method  uses  a 
quadratic  model  Hessian,  and  also  has  a  built-in  line  search  with  a  sufficient  decrease  condition.  The  line 
search  is  based  on  the  predicted  decrease  computed  from  the  quadratic  model.  If  the  actual  improvement 
of  the  objective  function,  J,  is  close  to  the  amount  predicted  by  the  model  Hessian  after  a  step  is  taken, 
then  the  method  decreases  the  Levenburg-Marquardt  parameter,  vc ,  effectively  increasing  the  relative  size 
of  the  next  step,  which  hopefully  accelerates  the  convergence.  As  vc  is  decreased  to  0  the  method  becomes 
Damped  Gauss-Newton  (meaning  Gauss-Newton  with  a  line  search).  If,  however,  the  actual  improvement 
of  J  after  a  step  is  not  sufficient  (or  is  even  negative),  vc  is  increased,  effectively  scaling  back  the  Newton 
step,  and  we  retest.  If  there  are  too  many  reductions  then  we  declare  a  “line  search  failure”  meaning  that 
too  small  a  step  is  required  to  decrease  the  objective  function. 

Usually  a  method  would  exit  after  a  line  search  failure,  returning  the  best  approximation  so  far.  But 
we  use  this  failure  to  call  an  adaptive  mesh  size  routine,  i.e. ,  an  Implicit  Filtering  technique.  The  idea  is 
that  the  failure  is  likely  due  to  the  fact  that  the  direction  the  finite  difference  gradient  chose  is  probably 
not  an  actual  “descent  direction”  in  the  global  sense.  In  other  words,  the  finite  differencing  is  most  likely 
differentiating  noise.  In  the  same  manner  that  a  smooth  surface  may  look  rough  under  a  microscope,  using 
too  small  of  a  differencing  step  amplifies  effects  from  round-off  error  and  other  sources  of  numerical  noise. 
Our  technique  is  to  increase  the  relative  differencing  step,  h,  recompute  the  gradients,  and  then  try  the 
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Levenburg-Marquardt  method  again.  The  relative  differencing  step,  h,  is  such  that  the  gradient,  V/t,  of 
J(q)  =  J([d,S))  is  computed  with 


V^([d,<S]) 


j((l  +  h)d,8)-J(d,S) 
hd 

j(d,(l+h)8)-j(d,8) 

hS 


We  apply  a  similar  approach  to  modifying  the  differencing  step  h  as  we  do  for  changing  vc  in  that  after  a 
successful  step  we  decrease  h,  but  if  we  have  another  failure  we  increase  h  even  more.  Since  the  convergence 
rates  of  gradient  based  methods  are  dependent  on  the  size  of  h  (for  example  Gauss-Newton  is  0(K2)),  we 
want  h  to  be  as  small  as  possible  and  still  be  effective,  similarly  with  uc.  We  use  a  three  tiered  approach 
to  changing  h.  Initially  we  set  h  =  10~9.  To  increase  h  we  raise  it  to  the  |  power,  to  decrease  we  raise  it 
to  the  |  power.  Additionally  we  define  ICG4  to  be  the  maximum  allowable  differencing  step  value.  Thus 
h  e  {ICG9,  ICG6,  ICG4}. 

In  general  an  optimization  method  exits  with  “success”  if  the  norm  of  the  current  gradient  is  less  than 
tol  times  the  norm  of  the  initial  gradient.  However,  in  our  method  we  do  not  immediately  trust  the  finite 
difference  gradients,  and  instead  call  Implicit  Filtering  when  the  gradients  appear  small.  When  we  have 
verified  small  gradients  on  all  three  scales  (the  various  values  of  the  differencing  step  h  defined  above),  then 
we  exit  with  “success”. 


Remark  4  In  practice,  a  very  good  solution  is  found  within  a  couple  of  Levenberg-Marquardt  steps,  and 
then  an  equal  number  of  Implicit  Filtering  iterations  verify,  and  sometimes  enhance,  this  solution.  In  the 
interest  of  efficiency,  and  since  this  is  a  parameter  identification  problem,  we  exit  early  with  “ success ”  if 
our  objective  function  is  satisfactorily  small  (i.e.,  tol  times  the  initial  value),  thus  saving  at  least  half  of  the 
possible  iterations. 

Additionally  we  impose  a  restriction  on  the  number  of  “pullbacks”  on  each  linesearch,  and  on  the  number 
of  linesearches,  effectively  limiting  the  total  number  of  iterations.  If  a  small  gradient  has  not  been  verified 
before  exhausting  the  maximum  number  of  iterations,  we  exit  with  “failure” . 


5.4  Numerical  Issues 

For  small  N  the  difficult  cases  are  those  with  large  depth.  This  is  because,  as  we  mentioned  before,  the 
computational  domain  is  effectively  increased  when  the  depth  is  increased,  making  the  mesh  sizes  larger  and 
increasing  the  level  of  numerical  error.  The  magnitude  of  <5  does  not  seem  to  have  a  significant  effect  on  the 
convergence  of  the  method. 

An  obvious  disadvantage  to  having  a  large  N  is  that  each  simulation  takes  much  longer.  Although  on 
average  fewer  function  calls  are  required  with  the  larger  N  cases  to  get  the  same  level  of  accuracy,  in  general 
the  total  execution  time  is  quadrupled  when  the  number  of  elements  is  doubled.  This  is  consistent  with  the 
fact  that  complexity  of  the  most  time  consuming  part  of  the  simulation,  the  linear  solves,  is  O(N),  and  the 
number  of  time  steps  Nt  is  also  O(N).  So  when  we  double  the  number  of  finite  elements  we  are  also  doubling 
the  number  of  time  steps.  Therefore,  we  get  an  overall  complexity  of  0(N2).  Thus,  as  mentioned  before,  in 
our  inverse  problem  we  choose  to  use  the  number  of  elements  just  above  the  threshold  of  when  numerical 
error  is  apparent. 

We  should  also  mention  that  in  order  to  create  data,  in  lieu  of  actual  experimental  data,  we  perform  a 
simulation  at  a  higher  resolution  believing  it  to  be  more  accurate.  Specifically,  we  double  the  number  of 
finite  elements.  Since  the  time  step,  and  therefore  the  effective  sample  rate  if  the  time  step  is  too  large,  are 
both  dependent  upon  the  mesh  size  (recall  At  =  O(h)  and  Remark  2),  the  sample  times  of  the  simulated 
data  do  not  necessarily  correspond  with  the  sample  times  of  the  simulations  at  the  lower  resolution.  (In 
general  we  have  twice  as  many  samples  from  the  higher  resolution.)  Thus  in  order  to  compute  the  modified 
least  squares  error  between  the  two  vectors,  we  perform  a  linear  interpolation  of  the  simulated  data  onto 
the  sample  times  at  the  lower  resolution.  See  Figure  55.  Note  that  in  the  usual  case  where  we  simply  have 
twice  as  many  sample  points  from  the  higher  resolution  simulation,  we  are  in  effect  discarding  sample  points 
rather  than  doing  a  true  interpolation. 
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Figure  55:  Plotted  are  the  actual  simulated  data  (N  =  2048),  the  interpolation  of  the  simulated  data  onto 
the  low  resolution  sample  times  (N  =  1024),  the  result  of  the  minimization  routine  (N  =  1024),  and  a  low 
resolution  (N  =  1024)  simulation  using  the  exact  values  of  d  and  6. 


For  comparison  we  compute  the  low  resolution  simulation  using  the  values  d*  and  5*  (note  that  this  is 
not  the  same  as  taking  the  high  resolution  simulation  and  interpolating  it  onto  the  low  resolution  time  steps, 
which  we  actually  use  as  our  observed  data).  In  every  case  that  we  have  tested,  </,  when  computed  with 
the  d  and  5  values  found  from  the  optimization  routine  ( dmin  and  6min ),  is  less  than  or  equal  to  J  when 
computed  with  the  original  values  (d*  and  5*).  This  suggests  that  an  actual  global  minimum  of  the  objective 
function  has  been  found,  even  though  the  final  estimates  of  d  and  6  themselves  are  not  necessarily  equal 
to  d*  and  6*.  Note  in  Figure  55  that  the  simulation  using  original  values,  (d*,S*),  is  in  fact  closer  to  the 
original  data,  but  the  simulation  using  the  minimizer  values,  (dm»n,5mjn),  is  closer  to  the  interpolated  data 
(see  for  example  the  [.335,  .3352]  interval). 

Although  we  could  compute  our  optimization  routine  at  the  same  resolution  as  the  simulated  data  to 
get  a  better  fit  in  our  tests,  this  would  not  properly  represent  the  real-life  phenomenon  of  sampling  data. 
Sampled  data  is  inherently  not  a  completely  accurate  representation  of  a  physical  observation.  We  believe 
that  our  interpolation  approach  gives  a  more  realistic  expectation  of  how  our  method  would  perform  given 
actual  experimental  data.  In  order  to  further  test  the  robustness  of  our  inverse  problem  solution  method  we 
introduce  random  noise  to  the  detected  data  in  Section  5.5. 

5.5  Numerical  Results 

Tables  4  and  5  show  the  final  computed  approximations  for  the  depth  of  the  slab  ( dmin )  and  the  width  of 
the  gap  behind  it  (Smin).  The  relative  differences  from  the  original  values  used  to  generate  the  data  (d*  and 
(5*),  are:  for  depth,  on  the  order  of  .0001  and  for  5,  on  the  order  of  .01.  However,  this  does  not  imply  that 
the  optimization  routine  was  unable  to  find  the  optimal  solution.  Recall  that  since  our  data  is  generated 
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with  essentially  a  different  simulator  than  our  toward  solves,  the  original  values  do  not  necessarily  minimize 
the  objective  function.  The  objective  function  values  give  a  better  indicator  of  how  well  the  optimization 
routine  works  since  it  shows  the  fit  to  the  generated  data.  Table  6  shows  the  final  objective  function  values. 
In  each  of  these  cases,  the  final  objective  function  value  ( Jmin )  was  less  than  J*  :=  J(q*).  In  fact,  the  ratios 
Jr  :=  Jmin/ J*  were  on  average  .3008.  We  consider  any  Jr  <  1  to  represent  a  successful  convergence. 

Although  6  values  that  are  near  |  =  3.7475  x  10_4m  are  the  most  difficult  for  which  to  obtain  initial 
approximations,  we  see  that  the  objective  function  values  in  these  cases  are  just  as  small  (and  the  final 
estimates  are  just  as  close)  as  for  other  S  values. 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

0.0200053 

0.0200022 

0.0200006 

0.0200005 

0.0200002 

.04 

(N=2048) 

0.0399948 

0.0399974 

0.0400005 

0.0400005 

0.0399999 

.08 

(N=4096) 

0.0799973 

0.0799987 

0.0800006 

0.0800006 

0.0800003 

.1 

(N=8192) 

0.0999945 

0.0999974 

0.1 

0.1 

0.0999999 

.2 

(N=16384) 

0.200011 

0.200005 

0.2 

0.2 

0.200001 

Table  4:  The  final  estimates  of  d. 
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d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.40622e-05 

0.000196754 

0.000398642 

0.000597275 

0.00079707 

.04 

(N=2048) 

0.000106435 

0.000203916 

0.000394204 

0.000592156 

0.000793622 

.08 

(N=4096) 

0.000103585 

0.000202273 

0.000395791 

0.000593861 

0.000794401 

.1 

(N=8192) 

0.000106593 

0.000203876 

0.000396203 

0.000594976 

0.000795985 

.2 

(N=16384) 

8.7456e-05 

0.000191808 

0.00040297 

0.000602902 

0.00080129 

Table  5:  The  final  estimates  of  6. 
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d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

0.00786171 

0.00906699 

0.0115657 

0.0233783 

0.0447687 

.04 

(N=2048) 

0.021516 

0.0343314 

0.0514108 

0.0700747 

0.0927117 

.08 

(N=4096) 

0.0116105 

0.0145428 

0.0201004 

0.0272513 

0.0344458 

.1 

(N=8192) 

0.00304723 

0.00547532 

0.00779186 

0.00931778 

0.0118529 

.2 

(N=  16384) 

0.000609258 

0.00133978 

0.00146975 

0.000962975 

0.000766141 

Table  6:  The  objective  function  value  of  the  final  estimates. 


The  execution  time,  in  seconds,  is  given  in  Table  7.  While  the  previous  tables  show  that  we  were  actually 
able  to  resolve  the  case  of  20cm  depth,  we  see  here  what  price  we  had  to  pay.  The  average  execution  times 
for  each  of  different  mesh  sizes  (N  =  1024,2048,4096,8192,  and  16384)  were  39,248,1452,6229,  and  35509 
seconds,  respectively.  Each  represents  and  increase  in  time  over  the  previous  mesh  size  by  a  factor  of 
6.4,5.9,4.28,  and  5.7,  respectively.  This  is  consistent  with  the  fact  that  the  forward  solves  are  order  0(h2). 
However,  the  additional  sample  points  for  the  larger  N  cases  allowed  for  smaller  objective  function  values 
which  took  increasingly  more  iterations  to  satisfy  the  relative  tolerance  in  our  stopping  criteria.  This  explains 
why  we  do  not  see  ratios  closer  to  the  expected  4  for  order  0(h2)  methods.  The  number  of  calls  to  the 
simulator  (related  to  the  number  of  Levenburg-Marquardt  iterations)  for  each  case  is  given  in  Table  8. 

5.5.1  Relative  Random  Noise 

We  add  random  noise  to  the  signal,  as  mentioned  above,  in  order  to  more  closely  simulate  the  experimental 
process  in  data  collection.  As  in  Section  4.3.2,  we  start  with  relative  noise  where  the  absolute  value  of  the 
noise  is  proportional  to  the  size  of  the  signal.  If  Et  is  the  data  sampled,  then  we  define  £)  =  E,(l  +  vrr]i), 
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d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

39 

35 

48 

38 

35 

.04 

(N=2048) 

198 

197 

319 

196 

328 

.08 

(N=4096) 

1560 

1458 

1500 

1071 

1671 

.1 

(N=8192) 

6141 

5256 

7625 

5319 

6802 

.2 

(N=16384) 

31818 

42755 

35860 

26521 

40593 

ible  7: 

The  execution  time  in 

seconds  to  find 

the  final  estimate 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

23 

23 

33 

26 

23 

.04 

(N=2048) 

24 

24 

43 

24 

45 

.08 

(N=4096) 

38 

35 

36 

24 

41 

.1 

(N=8192) 

29 

24 

37 

24 

33 

.2 

(N=16384) 

36 

50 

41 

29 

48 

Table  8:  The  total  number  of  calls  to  the  simulator  required  to  find  the  final  estimates. 


where  77 ,  are  independent  normally  distributed  random  variables  with  mean  zero  and  variance  one.  Again, 
the  coefficient  vr  determines  the  relative  magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  Ei,  in 
particular,  vr  =  0.01  corresponds  to  2%  noise.  We  tested  relative  magnitude  levels  of  2%,  10%,  and  20% 
(corresponding  to  vr  =  .01,. 05,  and  .1  respectively).  Table  9  shows  the  initial  estimates  of  d  for  vr  =  .1 
(other  vr  values  produced  very  simliar  results,  so  we  only  give  the  worst  case  here).  Tables  10  through  12 
show  the  initial  estimates  of  6  for  several  different  noise  levels.  There  is  not  a  noticable  significant  change  in 
the  accuracy  of  the  estimates  even  for  vr  =  .1.  In  nearly  all  the  cases  the  estimate  was  close  enough  for  the 
optimization  method  to  converge  (Jr  <  1)  to  the  expected  minimum.  The  only  exceptions  were  with  v  =  .1 
and  S  =  .0004,  which  are  understandably  the  most  difficult  cases. 

The  final  approximations  dmin  and  5min  in  the  presence  of  noise  are  given  in  Tables  13  through  18.  Some 
approximations  with  high  noise  may  appear  to  be  better  approximations  than  some  with  little  or  no  noise. 
For  example,  with  S*  =  .0001,  d*  =  .04,  the  vr  =  .1  approximations  are  an  order  of  magnitude  closer  to  the 
original  values  than  the  vr  =  0  approximations.  This  is  not  to  say  that  the  noise  helps  the  approximation 
method.  Rather,  it  is  for  the  same  reason  that,  for  example,  as  shown  in  Figure  55,  the  actual  parameter 
values  produced  a  signal  farther  away  (in  the  Least  Squares  sense)  from  the  generated  data  than  a  signal 
computed  with  the  approximated  parameter  values.  The  resulting  objective  function  values  give  a  better 
indication  of  the  accuracy  of  the  approximation  to  the  data.  A  comparison  of  Tables  6  and  19  clearly  show 
that  the  data  without  noise  is  more  accurately  matched  by  its  approximations  than  those  with  noise. 


48 


d 

1024 

N 

2048 

4096 

.02 

0.020019 

0.0199883 

0.0199998 

.04 

0.0399689 

0.0399766 

0.0399919 

.08 

0.0799225 

0.0800069 

0.0799992 

Table  9:  The  initial  estimates  of  d  with  vr  —  .1.  The  values  in  bold  denote  the  values  used  to  approximate 
5  in  Table  12.  (Initial  estimates  of  d  using  other  vr  values  were  very  simliar  and  therefore  are  omitted.) 

<5 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.43492e-05 

0.000127736 

0.000611592 

0.000647568 

0.000791472 

.04 

(N=2048) 

0.000115174 

0.000209026 

0.000611592 

0.000647568 

0.000755496 

.08 

(N=4096) 

3.20852e-05 

0.000194313 

0.000611592 

0.000647568 

0.000827448 

Table  10:  The  initial  estimates  of  <5  with 

vr  =  .01. 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

4.34312e-05 

0.000129248 

0.000575616 

0.000575616 

0.000755496 

.04 

(N=2048) 

0.000115175 

0.000211701 

0.000611592 

0.000755496 

0.000755496 

.08 

(N=4096) 

1.76356e-05 

0.000194415 

0.000647568 

0.000647568 

0.000683544 

Table  11:  The  initial  estimates  of  <5  with 

vr  =  .05. 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

0.000107928 

0.000172127 

0.000575616 

0.000323784 

0.000611592 

.04 

(N=2048) 

9.99018e-05 

0.000216157 

0.000575616 

0.000611592 

0.000647568 

.08 

(N=4096) 

4.76434e-05 

0.000194762 

0.000611592 

0.000647568 

0.000755496 

Table  12:  The  initial  estimates  of  8  with  vr  =  .1. 


49 


6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

.04 

.08 

(N=1024) 

(N=2048) 

(N=4096) 

0.0200057 

0.0399944 

0.079997 

0.0200021 

0.0399971 

0.0799986 

0.0200006 

0.0400005 

0.0800004 

0.0200008 

0.0400005 

0.0800005 

0.0200001 

0.0399998 

0.0800003 

Table  13:  The  final  estimates  of  d  using 

vr  =  .01. 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

.04 

.08 

(N=1024) 

(N=2048) 

(N=4096) 

0.0200061 

0.0399951 

0.0799957 

0.0200009 

0.0399982 

0.0799971 

0.0200016 

0.0399993 

0.0799999 

0.0200013 

0.0400006 

0.0800004 

0.0200002 

0.0399997 

0.0800001 

Table  14:  The  final  estimates  of  d  using 

vr  =  .05. 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

.04 

.08 

(N=1024) 

(N=2048) 

(N=4096) 

0.0200059 

0.0399993 

0.0799982 

0.020006 

0.0399976 

0.079997 

0.0200017 

0.0399991 

0.0800002 

0.0199998 

0.0400021 

0.0799986 

0.0199974 

0.0399982 

0.0800013 

Table  15:  The  final  estimates  of  d  using  ur  =  .1. 
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6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.34307e-05 

0.000196985 

0.000398454 

0.000596308 

0.000797028 

.04 

(N=2048) 

0.000106823 

0.000204701 

0.000394249 

0.0005922 

0.000793984 

.08 

(N=4096) 

0.000103996 

0.000202579 

0.000396222 

0.000593539 

0.000794799 

Table  16: 

The  final  estimates  of  S  using 

ur  =  .01. 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.40065e-05 

0.000198868 

0.000396695 

0.000597252 

0.000794729 

.04 

(N=2048) 

0.000105547 

0.000203208 

0.000400938 

0.000592641 

0.000791008 

.08 

(N=4096) 

0.000105656 

0.000203992 

0.000396485 

0.000594405 

0.000793268 

Table  17: 

The  final  estimates  of  S  using 

-3 

II 

o 

Or 

d 

.0001 

<5 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.34486e-05 

0.000188052 

0.000392001 

0.000594642 

0.00079958 

.04 

(N=2048) 

0.000100242 

0.000204073 

0.000401208 

0.000586155 

0.000798636 

.08 

(N=4096) 

0.000101771 

0.000205295 

0.000400441 

0.000604353 

0.000798047 

Table  18:  The  final  estimates  of  (5  using  vr  =  .1. 


6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

3.89043 

4.0585 

1.64814 

3.46279 

4.1071 

.04 

(N=2048) 

1.29514 

1.42229 

1.47275 

1.5261 

1.61906 

.08 

(N=4096) 

0.779184 

0.717672 

0.848017 

0.647817 

0.674641 

Table  19:  The  objective  function  value  of  the  final  estimates  using  vr  =  .1. 
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5.5.2  Standard  Error  Analysis 

In  an  actual  inverse  problem  using  data  collected  by  experiment,  one  needs  to  have  confidence  intervals 
on  all  parameter  estimates.  We  will  apply  standard  error  techniques  to  an  Ordinary  Least  Squares  (OLS) 
formulation  of  our  problem  to  obtain  confidence  intervals  on  our  estimates.  In  order  to  rewrite  our  objective 
function  in  an  OLS  formulation,  we  define  y{t\q )  =  \E(t,0;q)\  to  be  our  estimate  to  y  =  \E\,  which  is  the 
data  we  are  trying  to  fit  by  determining  q  =  ( d1  d).  Now  it  is  clear  that  our  objective  function  can  be  written 
in  the  standard  OLS  form 

1  Na 

J (<?)  =  y{U\i)  -&l  • 

i—  1 

For  simplicity  of  terminology,  in  this  section  alone,  we  will  refer  to  \Ei\  as  the  data  and  to  | E(ti,  0;  q)\  as  the 
simulations. 

With  the  relative  random  noise  described  above  we  do  not  have  constant  variance,  as  is  demonstrated  in 
Figures  56  and  57.  Here  we  have  plotted  the  residual  r.j  :=  \E(ti,  0;  qoLs)\  ~  |-E)|  against  time,  tj,  and  also 
against  \E(ti,  0;  qoLs)\-  As  one  would  expect  with  noise  that  is  relative  in  size  to  the  signal  value,  we  have 
a  pattern  in  Figure  56  that  follows  the  pattern  of  the  original  signal.  Figure  57  demonstrates  the  fan  shape 
associated  with  noise  that  is  dependent  upon  the  size  of  the  signal,  i.e.,  non  constant  variance. 

Since  constant  variance  is  most  conveniently  assumed  in  standard  error  analysis,  we  further  consider 
estimates  obtained  from  an  inverse  problem  applied  to  data  with  constant  variance  random  noise  added.  In 
particular,  the  data  we  now  consider  is  generated  by 

Ej  =  E(tj,  0, ;  q*)  +  f3vrr]j 


where 

rjj  ~  Af  (0, 1) 

and  the  constant  /3  is  a  scaling  factor  chosen  simply  so  that  the  noise  level,  vr,  will  somewhat  correspond  to 
the  parameter  vr  used  in  the  previous  section  on  relative  noise.  Specifically,  (3  =  maxj  Ej/10  ensures  that  J* 
in  the  constant  variance  cases  is  on  the  same  order  of  magnitude  as  those  in  the  relative  noise  cases  above 
for  all  choices  of  d  and  S  that  we  have  considered. 

The  variance  of  this  data  is 


a2  =  £§[/?2i#$  =  (32v2r£%2\  =  « 

where  £§  denotes  the  expectation.  Therefore,  we  do  have  constant  variance.  Note  further  the  resulting  lack 
of  patterns  in  Figures  58  and  59.  The  suspicious  looking  phenomenon  of  many  points  on  the  line  E  =  0  is 
simply  because  in  the  original  data  E  is  very  close  to  zero  most  of  the  time.  (Recall  Figures  45  and  46  for 
example.)  Figure  60  demonstrates  graphically  the  difference  between  relative  noise  and  constant  variance 
noise.  The  relative  noise  case  is  particularly  difficult  in  our  inverse  problem  since  most  of  our  initial  estimates 
are  based  on  accurately  determining  the  peak  locations,  yet  this  is  exactly  where  most  of  the  relative  noise 
is  concentrated. 

With  constant  variance,  and  further,  assuming  that  each  ly  is  identically  independently  (normally)  dis¬ 
tributed,  we  have  that  (see  [4])  in  the  limit  as  Ns  — >  oo 

qoLS  ~  A f2  (qo,  cro  [£T  {q0)£(qo)\  • 

Here  £(q)  =  ^p-(g)  which  is  an  Ns  x  2  matrix  since  q  =  (d,6)  and  \E\  is  evaluated  at  Ns  sample  times. 
Also,  the  scale  parameter  <Tq  is  approximately  given  by 


(Jn  = 


AT,  -  2 


Ne 

^(|E(ti,0;g0)|-|^|)' 


i= 1 


In  the  above  equations,  qo  denotes  the  theorectical  “true”  value  of  the  parameter  that  best  describes  the 
system  from  which  the  data  is  taken.  Note  that  in  this  case,  this  is  not  necessarily  the  same  as  q*  since  the 
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method  used  to  generate  the  data  is  different  from  the  forward  solve  simulator.  Therefore  qo  is  generally 
unknown  even  in  examples  with  simulated  data. 

As  demonstrated  in  the  previous  sections,  our  qoLS  is  often  a  better  minimizer  than  even  the  original  value 
of  q*,  therefore  we  will  approximate  qo  in  the  above  equations  by  qoLS-  hr  particular,  if  we  denote  the  covari¬ 
ance  matrix  as  Co  =  crfi  \£T  {qo)£{qo)\  ,  then  we  will  approximate  Co  by  C  =  ctqLS  [£' T{qoLs)£(dOLs)\  , 

where 

1  Ns  2 

aOLS  =  at  2  ^  (\E(tj,0;qoLs)\  —  |-^»|)  • 

s  i=  1 

We  compute  a qLS  by  multiplying  our  Jmin  by  an  appropriate  conversion  factor,  since  they  are  defined  in  a 
similar  manner.  However,  in  order  to  compute  the  partial  derivatives  with  respect  to  d  and  <5  in  £  we  employ 
forward  differencing,  which  requires  an  additional  forward  simulation  for  each  qj.  For  q  =  qoLS  we  have,  for 
example 

c  _d\E\  U  \E(ti,0;[q1,q2\)\-\E(ti,0;[(l- hd)q1,q2])\ 

til  —  “X — - - ; - 

oq\  hdq1 

and  similarly  for  each  £a.  In  our  computations  we  used  the  relative  differencing  factor  of  hd  =  1  x  10-4. 
One  could  also  use  a  sensitivity  equations  approach  (e.g.,  see  [1]  and  the  references  therein),  but  since  the 
variational  equations  are  quite  difficult  to  solve  for  this  example,  we  choose  instead  to  approximate  the 
partials  with  respect  to  q  directly  with  our  simulations. 

We  also  need  to  point  out  that  while  taking  the  absolute  value  of  a  function  limits  differentiability  at  a 
small  number  of  points,  the  derivative  does  exist  almost  everywhere.  The  absolute  value  function  does  not 
change  the  magnitude  of  the  derivative  where  it  exists,  which  is  what  we  need  to  compute  the  dot  product 
of  £  with  itself.  By  using  finite  differences  to  estimate  derivatives,  we  are  essentially  under-estimating  at 
the  discontinuities.  Under-estimating  a  few  points  out  of  thousands  is  not  going  to  significantly  change  our 
covariance  matrix.  (Alternatively,  one  could  have  defined  the  objective  function  by  squaring  the  signals 
instead  of  taking  absolute  values  to  avoid  this  problem.  In  this  research  we  were  interested  in  comparing  Ji 
and  J2  in  previous  sections  above  and  changing  the  scale  of  E  by  squaring  it  would  have  prevented  this.) 

With  £  calculated,  we  can  now  evaluate  C  =  ct‘qLS  \£t (qoLs)£(qoLs)\  1  ■  Then  the  standard  error  for 
qi  =  d  is  estimated  by  \/C\ 1  while  the  standard  error  for  q2  =  S  is  estimated  by  \JC22.  See  Tables  20  through 
27  for  confidence  intervals  relating  to  various  d*,  6*  and  vr  values.  For  example,  in  the  case  of  d*  =  .02, 
5*  =  .0002  and  with  vr  =  .01  our  covariance  matrix  is 


'  2.37122  x  10“15  -4.43815  x  10~15' 

-4.43815  x  10"15  9.1829  x  10"15 


which  results  in  the  confidence  intervals  d  €  (2.00004±4.86952  x  10  6)  x  10  2  and  6  £  (1.9941±0. 000958274)  x 
10"4. 

The  width  of  these  bounds  are  ±0.000243471%  and  ±0.0480555%  of  the  approximation  value  respectively. 
For  the  d*  =  .02  case,  the  average  size  of  the  confidence  intervals  for  vr  =  .01,  .05,  .1  respectively  were 
±.0002%,  ±.001%,  ±.002%  (averaged  over  various  <5*  values  ranging  from  .0001  to  .0008).  It  is  interesting  that 
the  widths  of  the  confidence  intervals  nearly  exactly  double,  on  average,  when  the  noise  level  doubles.  For  the 
d*  =  .04  case  the  average  size  of  the  confidence  intervals  were  ±.0001%,  ±.0006%,  ±.001%.  Likewise,  when 
the  widths  of  the  confidence  intervals  for  S*  =  .0002  are  averaged  over  several  various  d*  values  (.02,  .04,  .08) 
we  get  ±.05999%,  ±.2883%,  ±.5718%  for  vr  —  .01,  .05,  .1  respectively.  For  d*  =  .0004  the  averages  are 
±.03331%,  ±.1575%,  ±.3154%.  In  general,  larger  d*  and  8*  values  have  smaller  (tighter)  confidence  intervals. 
This  suggests  that  the  approximations  found  in  these  cases  are  better  than  those  estimating  small  parameters. 
While  this  is  intuitive,  it  is  not  apparent  looking  at  the  estimates  themselves  or  even  the  final  objective 
function  values  (see,  for  example,  Table  6). 
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|r(q)|vs  t,  (r(q)=  |e(q)|-|ehat|) 


Figure  56:  Plots  of  the  absolute  value  of  the  residual  r*  =  | E(ti,  0;  qoLs) I  —  \Ei\  versus  time  £,  when  the  data 
contains  relative  random  noise. 


|r(q)|  vs  |e(q)|,  (r(q)=  |e(q)|-]ehat|) 


Figure  57:  Plots  of  the  absolute  value  of  the  residual  rj  =  | E(ti,  0;  qoLs) I  —  \E,\  versus  the  absolute  value  of 
the  electric  field  E(ti,  0;  qoLS )  when  the  data  contains  relative  random  noise. 
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|r(q)|  vs  t,  (r(q)=  |e(q)|-|ehat|) 


t  (ns) 


Figure  58:  Plots  of  the  absolute  value  of  the  residual  =  \E(ti,  0;  qoLs) I  —  \Ei\  versus  time  ti  when  the  data 
contains  constant  variance  random  noise. 


|r(q)|  vs  |e(q)|,  (r(q)=  |e(q)|-|ehat|) 


Figure  59:  Plots  of  the  absolute  value  of  the  residual  r.;  =  \E(ti1 0;  qoLs) I  —  \Ei\  versus  the  absolute  value  of 
the  electric  field  E(ti,0;  qoLs)  when  the  data  contains  constant  variance  random  noise. 
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electric  field  (volts/meter) 


Relative  Noise  vs  Constant  Variance 


Figure  60:  The  difference  between  data  with  relative  noise  added  and  data  with  constant  variance  noise 
added  is  clearly  evident  when  E  is  close  to  zero  or  very  large. 
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6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00005  ±  9.30284  x  KT7)  x  10 
(2.00001  ±  6.50411  x  KT7)  x  10 
(2.00001  ±4.91232  x  10"7)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00013  ±  1.62162  x  10~6)  x  KT2 
(4.00001  ±  1.19064  x  10"6)  x  10"2 
(4.00002  ±  9.05240  x  10"7)  x  10"2 


Table  20:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  no  noise  (i.e., 
vr  =  0.0). 


6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00004  ±4.86952  x  KT6)  x  10 
(2.00001  ±  3.50259  x  10"6)  x  10 
(2.00001  ±  2.87772  x  10"6)  x  10 


~ 

-2 

-2 


d*  =  .04  ( N  =  4096) 

(4.00013  ±  5.69385  x  10~6)  x  KT2 
(4.00001  ±  4.02428  x  10"6)  x  10"2 
(4.00001  ±  3.32933  x  10"6)  x  KT2 


Table  21:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  vr  =  .01. 


6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00004  ±2.41541  x  10~5)  x  10 
(2.00000  ±  1.68896  x  KT5)  x  10 
(2.00003  ±  1.40398  x  10"5)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00014  ±  2.76640  x  10'5)  x  KT2 
(4.00001  ±  1.90853  x  10"5)  x  10"2 
(4.00000  ±  1.60390  x  10"5)  x  10"2 


Table  22:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  vr  =  .05. 


6 

.0002 

.0004 

.0008 


d* 

(2 

(2 

(2 


=  .02  (N  =  2048) 
.00000  ±4.72903  x 
.00003  ±  3.39327  x 
.00003  ±2.79911  x 


10"5)  x  10"2 
10"5)  x  10"2 
10"5)  x  10"2 


d*  =  .04  (N  =  4096) 

(4.00014  ±  5.48283  x  10~5)  x  10"2 
(4.00002  ±  3.87474  x  10"5)  x  10"2 
(4.00003  ±  3.19526  x  10"5)  x  10"2 


Table  23:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  vr  =  .1. 
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6  d*  =  .02  (N  =  2048)  d*  =  .04  (N  =  4096) 

.0002  (1.99272  ±  0.000182978)  x  KT4  (1.98142  ±  0.000317616)  x  10"4 

.0004  (4.00035  ±  0.000201885)  x  10"4  (4.00737  ±  0.000369841)  x  10"4 

.0008  (7.99833  ±  0.000136586)  x  10"4  (8.00332  ±  0.000251291)  x  10"4 

Table  24:  Confidence  intervals  for  the  OLS  estimate  of  6  when  the  data  is  generated  with  no  noise  (i.e., 
vr  =  0.0). 


<5 

d*  =  .02  (N  =  2048) 

d*  =  .04  (N  =  4096) 

.0002 

.0004 

.0008 

(1.99410  ±  0.000958274)  x  10"4 
(4.00170  ±  0.00108740)  x  10"4 
(7.99882  ±  0.000800042)  x  10“4 

(1.98029  ±  0.00111475)  x  10"4 
(4.00667  ±  0.0012499)  x  10"4 
(8.00486  ±  0.000923838)  x  10"4 

Table  25:  Confidence  intervals  for  the  OLS  estimate  of  <5  when  the  data  is  generated  with  noise  level  vr  =  .01. 


<5  d*  =  .02  (N  =  2048)  d*  =  .04  (N  =  4096) 

.0002  (1.99606  ±  0.00475672)  x  KT4  (1.98106  ±  0.00541764)  x  KT4 

.0004  (4.00190  ±  0.00524360)  x  10"4  (4.01214  ±  0.00593246)  x  10"4 

.0008  (7.99045  ±  0.00391181)  x  10"4  (8.00947  ±  0.00444525)  x  KT4 

Table  26:  Confidence  intervals  for  the  OLS  estimate  of  6  when  the  data  is  generated  with  noise  level  vr  =  .05. 


S 

d*  =  .02  (N  =  2048) 

d*  =  .04  (N  =  4096) 

.0002 

.0004 

.0008 

(2.00017  ±  0.00932701)  x  10“4 
(4.00070  ±  0.0105331)  x  10"4 
(7.99698  ±  0.00778563)  x  10"4 

(1.97674  ±  0.0107203)  x  10"4 
(4.01229  ±  0.0120445)  x  10"4 
(8.00361  ±  0.00886925)  x  10"4 

Table  27:  Confidence  intervals  for  the  OLS  estimate  of  <5  when  the  data  is  generated  with  noise  level  vr  =  .1. 
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6  Conclusion 


In  this  presentation,  we  have  explored  a  “proof  of  concept”  formulation  of  an  inverse  problem  to  detect  and 
characterize  voids  or  gaps  inside  of,  or  behind,  a  Dielectric  medium.  We  have  simplified  the  problem  to 
one  dimension  and  used  the  Maxwell  equations  to  model  a  pulsed  THz  normally  incident  electromagnetic 
interrogating  signal.  We  use  Finite  Element  discretization  in  space,  and  Finite  Differences  in  time,  to  simulate 
the  electric  field  in  the  time  domain.  This  is  coupled  with  a  Levenberg-Marquardt  scheme  in  an  optimization 
step  with  an  innovative  cost  functional  appropriate  for  reflected  waves  where  phase  differences  can  produce 
ill-posedness  in  the  inverse  problem  when  one  uses  the  usual  ordinary  least  squares  criterion.  We  have 
successfully  shown  that  it  is  possible  to  resolve  gap  widths  on  the  order  of  .2 mm  between  a  dielectric  slab 
of  20cm  and  a  metal  (perfectly  conducting)  surface. 

Future  work  on  this  problem  will  likely  involve  more  efficient  computational  methods  since  currently  the 
inverse  problem  involving  a  20cm.  slab  takes  10  hours.  Further,  more  sophisticated  models  for  describing 
the  polarization  mechanisms  in  non-homogenous  materials  must  be  developed.  Finally,  in  order  to  take 
scattering  and  non-normally  incident  eletromagnetic  signals  into  account,  multi-dimensional  models  will  be 
necessary. 
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